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ASPECTS  OF  THE  CALCULATION  OF 
BOUNDARY  LAYERS  ON  SHIPS'  PROPELLERS 
(FINAL  REPORT) 

EXECUTIVE  SUMMARY 

The  calculation  of  the  boundary  layer  on  marine 
propellers  is  important  for  several  reasons.  Viscous  losses 
in  propeller  performance  can  be  determined  by  such  a  computa¬ 
tional  capability.  However,  the  determination  of  the 
detailed  viscous  flow  in  the  propeller  root  and  tip  regions 
is  even  more  important  because  the  flow  in  these  regions 
directly  affects  the  hub  and  tip  vortex  generation  process. 
These  vortices  play  an  important  role  in  propeller  noise, 
cavitation  and  overall  performance.  Furthermore,  the  laminar 
boundary  layer  at  the  leading  edge  of  model  propellers  plays 
a  very  important  role  in  the  correlation  of  the  cavitation 
performance  of  model  to  full  scale  propellers. 

Computational  methods  for  calculating  the  boundary 
layer  on  the  faces  of  a  propeller  blade  have  been  developed 
by  Groves  (1981)  and  Groves  and  Chang  (1985).  These  methods 
ignore  or  only  crudely  approximate  the  flow  at  the  root  and 
leading  edge  of  the  propeller.  The  leading  edge  flow  has  a 
profound  effect  on  the  tip  flow  when  the  blades  are  swept 
back  as  on  highly  skewed  propellers. 


This  study  focussed  on  those  aspects  of  the  pro¬ 
peller  problem  that  have  not  been  treated  satisfactorily  by 
Groves  and  Chang.  In  particular  the  leading  edge  boundary 
layer  was  the  main  object  of  study.  However,  in  the  course 
of  this  study  another  important  aspect  of  leading  edge 
boundary  layers  emerged  that  has  considerably  wider  applica¬ 
tion  than  marine  propellers.  This  aspect  concerns  the  way 
boundary  layers  emerge  from  points  of  attachment  and  the  flow 
along  the  attachment  line  connecting  multiple  points  of 
attachment.  These  aspects  will  be  clarified  below.  At  pre¬ 
sent  it  need  only  be  noted  that  these  attachment  line  prob¬ 
lems  occur  also  at  the  bow  of  destroyer  hulls  with  a  pro¬ 
truding  sonar  bulb  and  also  at  wing  (or  foil)/body  inter¬ 
sections  if  the  wing  is  at  least  slightly  swept  foreward. 
Also,  the  flow  along  such  an  attachment  line  is  an  in¬ 
teresting  model  problem  that  may  have  important  implications 
to  the  horseshoe  vortex  rollup  at  the  leading  edge  of  a  wing 
body  intersection. 

This  report  is  divided  into  two  main  parts  de¬ 
scribing  the  research  performed  under  this  contract.  The 
first  part  describes  the  work  concerning  the  flow  on  the 
attachment  line  connecting  a  nodal  and  saddle  point  of  at¬ 
tachment.  The  second  part  describes  the  work  concerning  the 
flow  near  a  leading  edge  of  a  wing  or  blade. 
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The  research  performed  in  the  first  part, 


con¬ 


cerning  the  flow  between  nodal  and  saddle  points  of  at¬ 
tachment  had  a  very  interesting  outcome.  One  of  the  main 
findings  is  that  there  can  exist  separation  along  the  at¬ 
tachment  line  that  manifests  itself  as  a  horseshoe  vortex 
perpendicular  to  the  attachment  line  and  wraps  around  the 
body.  Under  certain  circumstances  this  horseshoe  vortex  is 
embedded  in  the  boundary  layer  and  the  separated  flow  does 
not  disturb  the  external  potential  flow.  Thus  this  flow  can 
be  calculated  by  ordinary  boundary  layer  theory  without 
viscous/inviscid  interactions.  The  various  aspects  of  this 
work  described  in  the  report  will  be  abridged  and  submitted 
to  technical  journals.  One  paper  has  already  been  submitted. 

The  second  part  of  this  research  concerning  the 
leading  edge  flow  on  blades  or  wings  has  not  been  as  suc¬ 
cessful.  A  model  problem  consisting  of  a  very  flat  oblate 
spheriod  (say  a  disk)  in  an  edgewise  stream  was  examined. 
The  boundary  layer  could  be  calculated  when  the  disk  was  at 
zero  angle  of  attack.  However,  great  difficulty  was  en¬ 
countered  in  developing  a  stable  method  that  would  properly 
account  for  the  flow  in  the  vicinity  of  the  curved  dividing 
streamline  that  separates  the  flow  going  over  the  top  and 
bottom  of  the  disk.  At  an  angle  of  attack  this  dividing 
streamline  is  no  longer  the  edge  symmetry  plane,  but  it  curls 
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off  the  edge  to  the  windward  (pressure)  side  of  the  disk.  A 
full  description  of  the  work  and  results  of  this  effort  is 
given  in  the  report.  We  give  below  a  summary  of  the  status 
of  the  propeller  leading  edge  problems 

(i)  A  numerical  method  for  dealing  with  the  3D  laminar 
boundary  layer  equations,  with  rotation,  and  in  non- 
orthogonal  surface  coordinates  has  been  developed  and 
programmed . 

(ii)  The  method  has  been  set  up  in  a  general  way  so  that  a 
variety  of  coordinate  sytems  can  be  handled  easily  as 
demonstrated  in  Sections  4  and  5  of  Part  II  of  this 
report . 

(iii)  The  method  can  deal  successfully  with  non  flat  bodies 
such  as  ship  bows  and  elongated  bodies  at  zero  or 
nonzero  angles  of  attack,  such  as  ovary  spheroids  or, 
more  practically,  submarines  or  torpedoes.  In  fact, 
the  same  basic  set  of  algorithms  and  codes  were  used 
for  the  direct  numerical  integration  of  the  attachment 
line  boundary  layer  of  Part  I.  This  type  of  calcula¬ 
tion  has  direct  relevance  to  destroyer  bows. 


(iv)  A  new  method  for  advancing  from  the  stagnation  point 
along  the  foreward  curve  of  attachment  near  a  leading 
edge  of  a  flattened  body  such  as  a  wing  or  blade  was 
unsuccessful.  The  method  does  not  seem  to  have  any 
flaws  and  at  this  time  it  is  not  certain  whether 
coding  errors  have  prevented  success  in  using  this 
method . 

(v)  Some  alternative  methods  to  achieve  success  in  cal¬ 
culating  the  leading  edge  flow  might  be  to  use  an 
explicit  backward  differencing  with  respect  to  the 
crossflow  coordinate  or  to  apply  the  characteristic 
box  method.  Resources  of  the  present  contract  have 
not  permitted  us  to  embark  on  efforts  to  develop  these 
alternative  methods.*  It  is  believed,  that  a  modest 
further  effort  can  resolve  the  problem  at  the 

attachment  line  and  result  in  a  very  useful,  and 
versatile  three-dimensional  boundary  layer  code. 


At  this  stage  only  our  inability  to  handle  the  starting 
conditions  at  a  general  attachment  line  prevents  our  code 
from  being  a  complete  general  3D  laminar  boundary  layer 
program. 
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C.  H.  von  Kerczek 


1. Introduction 

The  inviscid  surface  velocity  at  a  stagnation  point  of 
attachment  (x,y)=(0,0)=P  can  be  described  by  U»ax  and  V«*by  with  b>0 
in  the  orthogonal  coordinates  x  and  y  respectively.  The  stagnation 
point  of  attachment  P  is  characterized  by  the  value  of  c**a/b.  For  the 
values  of  c>  or  =0  the  point  P  is  a  nodal  point.  The  stagnation  point 
boundary  layer  at  a  nodal  point  was  completely  solved  by  Howarth 
(1951).  For  the  values  of  c<0  the  point  P  is  a  saddle  point.  The 
boundary  layer  at  a  saddle  point  was  solved  by  Davey(1961).  Davey 
showed  that  for  values  of  c<-1.0  no  boundary  layer  solutions  exist 
that  meet  the  outer  boundary  conditions  imposed  by  the  above  inviscid 
flow.  He  also  showed  that  for  values  of  c<-0.4294  the  flow  in  the  x 
direction  develops  a  reversed  flow  at  the  wall,  ie.,  the  boundary 
layer  separates.  However,  the  separated  flow  exists  as  a  solution  of 
the  stagnation  point  boundary  layer  equations  for  all  values  of  c  in 
the  range  -1.0<c<-0.4294. 

The  following  question  arizes:  Are  the  saddle  point  solutions 
affected  by  the  flow  along  the  symmetry  plane  (y*0)  that  emanates 
from  the  nodal  point?  Or,  is  the  saddle  point  flow  so  dominant  that 
no  matter  how  the  flow  begins  upstream  at  the  nodal  point,  the 
solution  always  evolves  into  Davey's  saddle  point  solution.  Cooke  and 
Robins  (1970)  investigated  this  question  by  numerically  integrating 


the  symmetry  plane  boundary  layer  equations  from  the  nodal  point  to 
the  saddle  point.  They  did  this  in  two  ways.  The  first  way  was  by 
series  expansions  about  the  nodal  and  saddle  point.  These  series  were 
found  to  overlap  and  thus  allowed  the  calculation  of  the  boundary 
layer  between  the  two  points  of  attachment.  The  second  way  they 
calculated  the  symmetry  plane  boundary  layer  was  by  direct  numerical 
integration. 

Cooke  and  Robins  carried  out  these  calculations  only  for 
values  of  c>or=-0.4294  so  that  the  boundary  layer  did  not  separate 
along  the  attachment  line.  Only  for  the  case  c»-0.4294  did  the 
boundary  layer  separate  at  the  saddle  point,  just  as  Davey  predicted 
it  would.  Thus  Cooke  and  Robins  showed  that  for  values  of  c>or- 
-0.4294  Davey's  saddle  point  solution  always  emerges  as  the  natural 
termination  of  a  flow  from  the  nodal  point  to  the  saddle  point.  The 
question  still  remains  if  the  separated  solutions  will  emerge  as  the 
natural  termination  of  the  flow  along  the  x-axis.  In  particular,  if 
the  flow  at  the  saddle  point  is  separated  then  there  must  be  a 
separation  point  on  the  x-axis  ahead  of  the  saddle  point.  Can  the 
boundary  layer  flow  proceed  through  this  separation  point  towards  the 
saddle  point  without  altering  the  external  inviscid  flow?  It  would 
seem  that  this  must  happen  because  Davey's  saddle  point  solutions 
exist  for  the  given  inviscid  flow.  Furthermore,  Klemp  and  Acrivos 
(1976)  have  shown  that  there  do  exist  boundary  layers  with  separation 
that  do  not  disturb  the  external  potential  flow. 

This  note  presents  some  calculations  that  show  the  attachment 
line  boundary  layer  between  a  nodal  and  saddle  point  for  some  values 
of  c<-0 . 4294  for  which  the  boundary  layer  along  the  x-axis  separates 
ahead  of  the  saddle  point.  Nevertheless,  the  saddle  point  is  reached 
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with  Davey's  separated  velocity  profile.  It  should  be  mentioned  that 
this  attachment  line  boundary  layer  finds  practical  application  to 
the  flow  along  the  stem  of  bulbous  bowed  ships. 

2. Formulation  of  the  Attachment  Line  Flow 

Cooke  and  Robins  considered  the  model  flow  over  a  plane 
surface  z-0  on  which  the  inviscid  velocity  (U,V,W)  is  given  by 

U»ax(l-x)+az2,  V=by,  W=-(a+b-2ax)z.  (1) 

The  z-axis  is  perpendicular  to  the  surface.  The  nodal  point  of 
attachment  is  the  point  (0,0,0)  and  the  saddle  point  is  (1,0,0).  The 
line  of  attachment  is  y=*0,  z=0.  The  boundary  layer  flow  velocity 

(u',v',w')  in  the  (x,y,z)  directions  respectively  is  rescaled  as 
follows: 

u'=Uu(x,z'),  v'=Vv(x,z'),  w’=(b/Re) 1^2w(x,z' )  (2) 

where  z *=(bRe) li/2z  and  Re  is  the  Reynolds  number.  The  resulting  y-0 
symmetry  plane  boundary  layer  equations  are 

w  ,+ex(l-x)u  +e(l-2x)u+v=0,  (3) 

Z  X 

u„,_,-wu_,-e(l-2x)u2-ex(l-x)uu  =-e(l-2x) ,  (4) 

z  Z  Z  X 

Vz*-wV-v2-ex(1-x)uV-1-  f5) 


Equation  (3)  is  the  continuity  equation  and  (4)  and  (5)  are 


the  x  and  y  momentum  equations  respectively.  The  parameter  e-a/b, 
b>0.  For  e>0  the  nodal  point  is  at  x-0  and  the  saddle  point  is  at 
x«l.  Thus  at  the  saddle  point  x»l  e  is  the  same  as  -c.  For  negative 
values  of  e  x*0  is  the  saddle  point  and  x=l  is  the  nodal  point.  For 
e>0  the  flow  along  the  x-axis  is  from  x«0  to  x=l  and  the  v-component 
of  the  flow  is  always  away  from  the  x-axis.  The  equations  (3-5)  are 
quasi  two-dimensional  and  are  solved  subject  to  the  boundary 
conditions 

u»v=w=0  at  z'*0  and  u=v=l  at  z'=infinity.  (6) 

These  equations  were  solved  numerically  using  the  well  known 
Keller  box  method  (see  Cebeci  and  Bradshaw, 1977) .  In  the  region  of 
flow  separation  the  numerical  integration  was  continued  in  the 
forward  direction  (x>0)  by  employing  the  FLARE  approximation.  This 
approximation  sets  u=0  whenever  u<0  and  u  occurs  as  a  coefficient  of 

V 

3. Results  and  Discussion 

The  results  of  the  numerical  solution  of  equations  (3-6)  are 

shown  in  Figure  1.  There  the  wall  skin  friction  coefficient  Cfx- 
1/2 

(bRe)  '  u(du/dz')0  is  plotted  versus  x  for  three  values  of  e.  The 
case  e-0.25  corresponds  to  no  separation  of  the  flow  in  the  x- 
direction.  This  case  was  presented  by  Cooke  and  Robins  and  the 

present  results  are  virtually  identical  to  theirs.  In  the  second  and 
third  cases  of  e»0.5  and  0.6  the  u  component  of  the  boundary  layer 
separates  ahead  (x<l)  of  the  saddle  point  of  attachment.  However  the 

numerical  integration  method  was  able  to  proceed  through  the 

separated  flow  region  to  the  saddle  point.  For  values  of  e  slightly 
larger  than  0.6  the  marching  method  breaks  down  before  the  saddle 
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point  is  reached.  The  values  of  Cfx/U  predicted  by  the  numerical 
integration  at  the  saddle  point  x«l  are  -0.1115  and  -0.2668  for  e»0.5 
and  0.6  respectively.  These  values  of  Cfx/U  are  in  very  good 
agreement  with  Davey's  results.  They  show  that  even  in  the  case  of 
separated  flow  Davey's  saddle  point  solutions  emerge  from  the  flow 
towards  the  saddle  point/  albeit  for  values  of  e  only  slightly  larger 
than  0.4294. 

The  question  remains  wether  flow  separation  along  the 
attachment  line  y«0  can  occur  without  disrupting  the  external 
inviscid  flow  for  larger  values  of  e.  Although  the  present  numerical 
integration  method  using  the  FLARE  approximation  failed  to  yield 
solutions  all  the  way  to  the  saddle  point  for  values  of  e>0.6  this  is 
not  conclusive  evidence  that  such  solutions  with  the  given  inviscid 
flow  are  not  possible.  Davey  speculates  that  once  separation  does 
occur  on  the  attachment  line  the  boundary  layer  flow  interacts  with 
the  external  inviscid  flow  in  such  a  way  that  the  saddle  point 
becomes  a  nodal  point  of  attachment.  This  certainly  must  occur  for 
e>l/  but  not  necessarily  for  e<l.  In  order  to  verify  if  separated 
solutions  with  the  inviscid  flow  (1)  are  possible  for  values  of  e>0.6 
a  more  robust  method  for  solving  these  boundary  layer  equations  is 
needed.  The  method  of  Klemp  and  Acrivos  (1976)  seems  to  be  a 
possibility.  They  treated  a  problem  in  which  the  boundary  layer 
containes  an  imbedded  separated  region  but  the  external  inviscid  flow 
was  unaffected  by  this  separation.  Another  method  that  is  presently 
being  persued  is  the  use  of  high  order  series  expansions  in  terms  of 
the  parameter  e.  The  point  of  expansion  is  the  two-dimensional 
stagnation  line  solution  given  by  e«0. 
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Figure  Caption: 

Fig.l  Plot  of  Cfx  vs  x  for  values  of  e-0.25,0.5  and  0.6.  The  nodal 
point  is  at  x«0  and  the  saddle  point  is  at  x-1. 
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Section  1 
INTRODUCTION 


Consider  the  viscous  flow  around  a  corrugated 
circular  cylinder  as  shown  in  Figure  1.  The  approaching  flow 
is  perpendicular  to  the  axis  of  the  cylinder.  Let  the  curve 
marked  by  the  points  labeled  N  and  S  be  the  intersection  of 
the  cylinder  with  the  meridian  plane  containing  the  cylinder 
axis  and  the  oncoming  stream  velocity  vector.  Let  e  be  a 
measure  of  the  amplitude  of  the  corrugations. 


The  curve  containing  the  points  N  and  S  is  the 
attachment  line  of  the  flow  and  the  points  N  and  S  are  the 
nodal  and  saddle  points  of  attachment  respectively.  The 
points  N  and  s  are  stagnation  points  of  the  external  in- 
viscid  flow.  At  the  attachment  line  the  flow  divides  and 
proceeds  around  the  circumference  of  the  cylinder.  In  the 
case  of  a  straight  cylinder,  e=0,  the  attachment  line  is 
entirely  one  of  stagnation.  It  is  then  the  classical  2D 
stagnation  line  found  in  all  planar  flows  about  cylindrical 
bodies.  In  the  case  of  the  corrugated  cylinder  there  is  a 
spanwise  flow  along  the  attachment  line.  This  spanwise  flow 
emanates  from  the  nodal  point  and  terminates  at  the  saddle 
point.  A  three  dimensional  boundary  layer  calculation  around 
the  corrugated  cylinder  must  begin  with  the  calculation  of 
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the  attachment  line  flow.  This  is  analogous  to  the  calcula¬ 
tion  of  the  symmetry  plane  boundary  layer  in  ordinary  three 
dimensional  boundary  layers. 

This  corrrugated  cylinder  boundary  layer  is  being 
considered  because  it  is  a  model  problem  for  several  im¬ 
portant  practical  boundary  layer  problems.  The  most  im¬ 
portant  of  these  practical  problems  is  the  flow  at  the  bow  of 
a  sonar  domed  destroyer  hull  as  shown  schematically  in  Figure 
2.  Boundary  layer  calculation  methods  must  start  at  stagna¬ 
tion  points  and  lines  of  attachment.  Up  to  the  present  time 
all  calculation  methods  have  only  dealt  with  simple  nodal 
points  of  attachment  in  general  three  dimensional  flows. 
Thus  these  methods  cannot  deal  properly  with  flow  around 
bulbous  bowed  destroyers.  Indeed,  practical  ship  boundary 
layer  calculation  methods  deal  with  turbulent  flows  that 
begin  downstream  of  the  front  stagnation  point.  The  calcula¬ 
tion  methods  are  started  with  artifical  initial  data  (or 
experimental  data  when  available,  see  Groves,  1981).  How¬ 
ever,  on  bow  geometries  as  the  one  depicted  in  Figure  2  the 
boundary  layer  emanating  from  the  attachment  line  between  the 
nodal  and  saddle  point  can  separate  to  form  a  horseshoe 
vortex  wrapped  around  the  stem  of  the  ship  and  trailing 
downstream  along  the  sides.  Thus,  ignoring  the  bow  flow 
would  result  in  missing  this  very  important  aspect  of  the 
boundary  layer. 


It  is  also  fairly  obvious  that  if  the  parameter  e 
characterizing  the  corrugation  amplitude  becomes  large  the 
geometry  will  have  a  certain  crude  resemblance  to  a  filleted 
hull/appendage  intersection.  This  type  of  intersection  prob¬ 
lem  is  very  important  in  many  flow  configurations.  In  the 
present  study  it  was,  in  fact,  the  hub/blade  intersection 
that  led  us  to  consider  the  corrugated  cylinder  problem. 

In  this  study  only  the  flow  along  the  attachment 
line  between  points  N  and  S  was  considered.  This  flow 
serves  as  the  starting  point  of  a  complete  calculation  of  the 
flow  around  the  cylinder.  In  the  overall  flow  field  the 
attachment  line  flow  serves  as  the  initial  conditions  for  the 
cross-flow.  The  mainstream  flow  is  around  the  circumference 
perpendicular  to  the  attachment  line. 

The  external  inviscid  flow  driving  the  attachment 
line  flow  and  the  governing  laminar  boundary  layer  equations 
are  given  in  Section  2.  Section  3  describes  the  flow 
directly  at  the  saddle  point  of  attachment.  This  flow  can  be 
considered  in  isolation.  Indeed,  Davey  (1961)  solved  the 
saddle  point  flow  completely,  but  left  many  unanswered 
questions  concerning  the  relationship  of  the  saddle  point 
flow  to  the  rest  of  the  attachmment  line  flow.  This  flow  has 
been  reconsidered  here  using  a  different  kind  of  numerical 


method  in  order  to  set  the  stage  for  the  calculation  of  the 
attachment  line  flow.  Section  4  details  the  results  and 
limitations  of  a  straightforward  numerical  integration  of  the 
attachment  line  flow.  Section  5  presents  a  high  order  per¬ 
turbation  method  to  calculate  the  attachment  line  flow  and 
the  results  obtained  by  this  method.  The  last  section. 
Section  6,  gives  some  closing  remarks  concerning  methods  for 
solving  the  attachment  line  flow  on  practical  configurations 
and  extension  to  turbulent  flows. 


Section  2 


FORMULATION  OF  THE  ATTACHMENT  LINE  FLOW 


Let  the  surface  arclength  coordinate  along  the 
attachment  line  be  denoted  by  x  and  let  y  be  the  coordinate 
normal  to  x  and  tangent  to  the  surface  of  the  cylinder.  Let 
z  be  the  coordinate  normal  to  the  cylinder  surface.  Then  in 
the  vicinity  of  the  attachment  line  the  potential  flow  on  the 

•f 

surface  of  the  cylinder  has  the  velocity  distribution  V  given 
by 


V  *  e  sinx  1  +  y  j  .  (1) 

The  surface  velocity  field  V  given  by  equation  (1) 
is  obtained  by  considering  small  amplitude  corrugations  and 
calculating  the  small  deviations  of  the  resulting  flow  from 
the  flow  around  a  straight  circular  cylinder  (see  Davey, 
1961).  There  are  an  infinite  number  of  nodal  points  at 
x=2nnx  and  an  infinite  number  of  saddle  points  at 
x=(2n+l)Ttx/2  for  n=0 ,  +1 ,  +2 , . . .  Only  the  node  at  x=0  and  the 
saddle  at  x=tt  is  considered  in  this  study.  The  rest  of  the 
nodes  and  saddles  are  the  same  as  x=0  and  tt  respectively  by 
symmetry. 


The  three  dimensional  boundary  layer  equations  in 
general  orthogonal  coordinates  are  the  following  (see 
Rosenhead,  1963):  The  continuity  equation  is 


—  +  ~  +  ~  +  K..U  +  K  v  =  0 
ax  ay  az  i  2 


The  momentum  equation  along  the  attachment  line  x  is 


u—  +  +  w—  +  K,  uv  -  K  v2  = 

ax  ay  az  i  2 


-  3E  +  (3) 

9x  3z2 


The  momentum  equation  in  the  y  direction  is 


u~  +  +  K  uv  -  k  u2  = 

3x  ay  az  2  1 


ap  ,  3  v 

-  +  V - =r 


The  boundary  layer  velocity  is  denoted  by  (u,v,w) 
in  the  ( x,y,z )  directions  respectively.  The  pressure  is 
denoted  by  p  and  can  be  computed  by  Bernoulli's  equation  as 


P0-p=(l/2)pV  •  V 


It  is  assumed  here  that  x,  y  and  z  are  arclength  coordinates 
so  that  the  metric  coefficients  can  be  suppressed. 


The  quantities  and  K 2  are  the  geodesic  curva¬ 
tures  of  the  x=constant  and  y=constant  curves  respectively. 
The  attachment  line  is  a  symmetry  line  so  the  K2=0.  The 
geodesic  curvature  K1={ 1/r )dr/dx,  where  r=r(x)  is  the  local 
radius  of  the  cylinder.  However,  is  neglected  in  this 
study  simply  because  it  plays  a  negligible  qualitative  role. 
The  term  can  be  easily  included,  but  in  this  model  problem 
only  the  main  pressure  distribution  is  of  interest. 


The  boundary  conditions  on  u,  v  and  w  are  the 


following : 


u  =  v  =  w  =  0 


at  z=0  , 


u  =  e  sin  x  ,  v»y  at  z  +°°  . 


The  boundary  layer  velocity  field  (u,v,w)  is  re¬ 


scaled  as 


u  =  e  sinxf  ,  v  =  yh  ,  w  =  Vv  g 


and  the  vertical  coordinate  z'  is  defined  as 


z  =>/v  z '  . 


Furthermore,  the  velocity  gradients  q  and  r  are  defined  as 


W 


VVMrJ; J 


and 


f  =  h  =  1  at  z  +  ”  .  (18) 

These  boundary  layer  equations  (12  -  18)  can  be 
solved  by  marching  from  the  nodal  point  x=*0  to  the  saddle 
point  x=tt.  Note  that  the  inviscid  flow  vanishes  at  the 
saddle  point  and  that  the  u  component  of  velocity  points  into 
the  saddle  point,  much  like  a  rear  stagnation  point.  It  is 
well  known  (see  Rosenhead,  1963)  that  the  boundary  layer 
equations  have  no  solutions  at  a  nodal  point  of  separation. 
But  a  saddle  point  differs  from  a  nodal  point  of  separation 
in  that  the  velocity  in  one  direction  flows  away  from  the 
stagnation  point.  This  outflow  is  so  strong  that  the  verti¬ 
cal  flow  at  the  stagnation  point  is  downward  rather  than 
upward.  For  example,  the  Inviscid  flow  Eq.  (1)  has  the  form 

V  =  e  sin  x  i  +  y  j  -  z  (1+ecosx)  k  (19) 

for  small  values  of  z  just  above  the  surface.  At  the  saddle 
point  x=ff  the  k  component  of  velocity  is  towards  the  body  as 
long  as  e<l  and  thus  the  boundary  layer  vorticity  is  not 
transported  away  from  the  wall.  Once  e>l  the  k  component  of 
velocity  is  away  from  the  wall  and  the  saddle  point  behaves 
like  a  nodal  point  of  separation  (i.e.,  a  rear  stagnation 
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point  on  a  closed  body)  at  which  steady  boundary  layer  solu¬ 
tions  do  not  exist. 


It  is  also  interesting  to  note  that  the  parameter  e 
characterizes  the  local  geometry  at  the  nodal  and  saddle 
points.  Por  example,  at  the  nodal  point  x=0,  e=l  corresponds 
to  a  spherical  stagnation  point  for  which  the  curvature  is 
the  same  in  every  direction.  The  value  e=0  is  the  two  dimen¬ 
sional  stagnation  line.  This  knowledge  suggests  an  alterna¬ 
tive  way  of  dealing  with  equations  (12  -  18).  This  alterna¬ 
tive  is  to  calculate  the  boundary  layer  as  an  asymptotic 
expansion  in  e  about  e=0.  This  is  indeed  one  of  the  tech¬ 
niques  used  in  the  present  study  and  will  be  described  in  a 
later  section. 

As  a  preliminary  analysis  of  the  full  attachment 
line  problem,  Eqs.  (12  -  18),  the  special  case  of  the  saddle 
point  flow  x=  it  will  be  examined  by  a  perturbation  expansion 
about  e=0  in  the  next  section. 


Section  3 


THE  SADDLE  POINT  FLOW 


3.1  BASIC  EQUATIONS 


The  flow  at  the  saddle  point  of  attachment,  x*»n  in 
Eqs .  (12  -  18),  is  governed  by  the  following  eqs. 


ah 

az 


+ 


r 


0 


(20) 


I3  +  h  =  ef 
az 


(21) 


3r 

3z 


gr  +  h 


2 


1 


3f  .  n 
s?  +  q  *  0 


*3  . 
az 


gq 


(22) 


(23) 


(24) 


with  the  boundary  conditi*  as 


f  *  g  »  h  »  0 


at  z=0 


(25) 


and 


f  =  h  =  1 


at  z=°°  . 


(26) 


In  fact  these  equations  also  govern  the  flow  at  the  nodal 
point  for  values  of  -l<e<0.  Davey  (1961)  solved  these  equa¬ 
tions  numerically  for  values  of  0<e<l.  He  found  that  for 
values  of  e>0.4294  the  flow  at  the  saddle  point  shows  a 
separated  u  (f)  component  velocity  profile.  This  indicates 
that  the  attachment  line  flow  will  separate  ahead  of  the 
saddle  point  for  values  of  e  in  the  range  0.4294<e<l.  This 
separation  is  what  initiates  the  horseshoe  vortex  that  will 
wrap  around  the  necked  down  region  of  the  corrugated 
cylinder.  However,  Davey  does  not  make  any  calculations  on 
the  attachmment  line. 


In  the  present  study  the  attachment  line  equations 
will  be  solved  by  direct  numerical  integration  and  by  small 
e-expansions.  For  this  reason  it  is  worthwhile  to  consider 
the  simpler  case  of  the  small  e-expansion  at  the  saddle  point 
to  check  the  results  with  Davey's  solutions.  Furthermore, 
this  small  e-expansion  at  the  saddle  point  will  indicate  the 
range  of  validity,  for  a  given  number  of  terms,  of  the 
general  case. 
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The  vector  (h,g,r,f,q) 


F  is  expanded  in  e  as 


follows: 


F 


CO 


2 


n 


F 

n 


where 


n 


=  F  (z)  =  (h, 
n  n 


V 


n 


fn' 


V 


By  substituting  Fn  into  equations  (20  -  26)  and 
terms  according  to  their  power  of  e,  the  following 
of  perturbation  equations  result. 


0(1): 


3h 

+  r  =  0 
0  Z  O 


3r 
_ o 

3z 


h  =  0 
o 


g  r  +  h  2 
o  o  o 


1 


> 


j 


(27) 


grouping 

sequence 


( 28a, bf c) 


o 


(28d,e) 


with  boundary  conditions 


0(e): 


*  0  at  z=0 

o  o  o 


ho  =  1 


at  z* 


3h 

— i  +  r  =0 
9z  1 

3g, 

—k —  +  h  —  f 
3  z  1  o 


3z 


-  <Vl  +  rogl  ■ 


2hohl) 


=  0 


(28f,g) 


( 29a,b,c ) 


0 


aq-, 

TT  '  'Vl  +  Vl1 


(29d,e) 


with  boundary  conditions 


at  z=0 

( 29f ,g) 

at  z*  . 


h  m  hi  =  9i  ■  0 

f.  *  h.  =  0 


in  the  general  case  of  0(e)  for  n>2  the  equations  are 


3h 

— 2-  +  r  =  0 
3z  n 


3g 

n 

3z 


+  h  = 
n 


'n-1 


>  (30a, b,c) 


3r 


n-1 


n 


Tz 


r  -  (g 


r  +  r  g  -  2h  h  ) 
o  n  o^n  o  n 


=  2 


(g  .  r 
yD  n-D 


-  h  .h  .) 
D  n-j 


j=l 


3  f 


a  +  q  =  0 


3  z  n 


3q 


n 


3z 


"  -  gQqn  =  J  g jqn_ j 

j-1 


with  boundary  conditions 


(30d,e) 


f  =  h  =  g„  =  0  at  z=o 

n  n  n 

f  =  h  =  0  at  z=°°  . 

n  n 


(30f,g) 


At  order  1  the  equations  (28a, b,c)  are  simply  the 
planar  stagnation  point  equations.  It  is  interesting  to  note 
that  in  this  case  the  equations  (28d,e)  are  exactly  the  same 


as  the  temperature  equation  with  Prandtl  number  equal  to  one 


at  a  planar  stagnation  point.  If  the  function  gQ(z)  is  known 

then  attachment  line  flow  (f^,q,J  (and  all  the  f  ,q  's  for 

o  o  n  n 

that  matter)  can  be  found  by  quadratures.  This  will  be  very 
useful  in  the  subsequent  analysis. 

3.2  NUMERICAL  METHOD 

The  problems  (28)  -  (30)  can  be  solved  sequentially 
by  numerical  methods.  It  is  desirable  to  calculate  as  many 
terms  as  possible  of  this  perturbation  series.  If  a  suf¬ 
ficient  number  of  terms  of  sufficiently  high  accuracy  can  be 
obtained  then  it  may  be  possible,  even  for  a  limited  radius 
of  convergence,  to  obtain  values  of  the  flow  velocities  over 
the  entire  range  of  values  of  e.  Thus,  the  numerical  method 
for  solving  the  differential  equations  in  problems  (28) 
(30)  needs  to  be  as  accurate  and  efficient  as  possible.  The 
following  numerical  method,  based  on  Hermite  interpolation, 
is  used.  This  method  will  be  called  the  "Hermite  box" 
method. 


The  equations  (28a, b,c)  -  (30a,,b,c)  all  reduce  to 

the  form 


~  =  F '  =  A  *  F  +  b 
dz 


. »  .  -  h  • 
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where  the  n  subscript  on  F  has  been  suppressed.  The  matrix  A 
can  be  a  function  of  F  and  thus  the  equations  are  nonlinear. 
However,  these  equations  can  be  quasi-linearized  in  a  simple 
form  such  as  replacing  A(F)  by  A(F*)  where  F*  is  the  value  of 
?  at  a  prior  iteration  step.  Then  an  iterative  procedure  of 
solving  (31)  with  A(F *)  will  converge  to  the  desired  solution 
F.  It  should  be  noted  that  only  the  equations  (28a, b,c)  are 


nonlinear.  Once  these  have  been  solved  then  the  rest  of  the 
equations  are  linear  and  the  coefficient  matrices  only  depend 


on  Fq.  Thus  only  the  equations  (28a, b,c)  need  to  be  quasi- 
linearized  and  iterated  to  convergence.  It  is  presumed  from 
here  on  that  this  is  done. 


By  Hermite  interpolation  it  can  be  shown  that  equa¬ 
tions  (31)  can  be  numerically  solved  on  a  partition  {z^  of 
the  z  interval  where  zro  >>  1  but  zm  <  °°,  by  any  of  the 


following  sequence  of  finite  difference  methods: 


+  Azi  -  Azi  +  , ,  ,rx 

F.  =  F.  .  +  (F.  '  +  F.  )  +  F"  (£) 

i  i-l  2  l  l-l  xz 


Az . 

F,  +  — —  <?  ' 
l-l  2  l 


+  +  F.^ 


Azi  -*(5) 

”>  +  720  F  W 


%»•  s*  V>*. 


Ill 


F. 

l 


(V 


(-1.  '  '  +  F.  ") 
l  l-l 


Azf 

l 

120 


Az ! 


(Fi*  '  ' 


T?  »  »  I 

Fi-1 


)  + 


100800 


?(7> 


(?)  (34) 


etc.  for  even  higher  orders.  In  these  formulas  primes  and 
superscripts  in  parenthesis  indicate  differentiation  with 
respect  to  z. 


The  derivatives  of  F  can  be  obtained  exactly  by 
using  the  original  differential  equation  (32).  For  example, 

F"  =  (A-F+S)’  *  A’  -  F  +  A  *  (A'P+b)  +  b* 


so  that 


F' '  =  (A2+A' )  *  F  +  b'  +  b  . 

The  higher  derivatives  are  obtained  in  an  analogous  manner. 
The  resulting  system  of  algebraic  equations  depend  only  on 
powers  and  derivatives  of  the  matrix  A  and  the  right  hand 
side  vector  b,  both  of  which  are  known  at  each  step  of  the 
calculation. 


Method  (32)  is  simply  the  standard  trapezoidal  rale 
of  integration  (now  referred  to  as  the  box  method).  Formulas 
(33)  and  (34)  are  higher  order  generalizations  of  the  trape- 


zoidal 

rule . 

Method 

(33)  is  now  often 

referred  to  as 

the 

"super 

box* 

method. 

We  prefer  to  call 

these 

Hermite 

box 

methods 

in 

general  and  "nth  order  Hermite  box" 

method 

for 

formulas 

using  the  nth 

derivative  of  F. 

In 

this  work 

formula  (34)  is 

used. 

It  has 

the 
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incredibly  small  truncation  error  (Az^  /100800)  F  (£)  where 
£  is  some  mean  value  constant.  All  these  formulas  can  be 
proved  to  be  absolutely  stable  so  that  large  step  sizes  Az^ 
can  be  fully  utilized  and  still  achieve  accurate  results.  In 
fact,  solution  of  equations  (28a, b,c)  for  the  planar  stagna¬ 
tion  point  produced  8-significant  figure  accuracy  with  50 
points  across  the  interval  z  belonging  to  [0,10].  This 
compares  favorably  to  Rosenhead’s  (1963)  solution  to  6  figure 
accuracy  which  required  about  50  points  in  the  interval  [ 0 , 5 J 
(about  half  our  step  size).  Rosenhead's  solution  used  a 
fourth  order  method  with  extrapolation  to  sixth  order. 

All  the  calculations  performed  in  this  work  used 

double  precision  arithmatic  and  formula  (34)  with  a  uniform 

step  size  of  Az  =  0.1.  Thus,  if  the  maximum  value  of 
-►(  7 ) 

|F  (£ )  |  is  0(1)  then  the  numerical  integration  errors  are 


about  10  .  Extensive  testing  on  various  step  sizes  and  at 

single  and  double  precision  confirms  this  estimate. 


CALCULATION  OF  THE  COMPONENTS  f  AND  q 


For  each  value  of  n  equations  (28d,e)  -  (30d,e)  all 
have  the  form 


3f 

— H  +  q  =  0 

3z 


(35a, b) 


TT  '  goqn  =  Sn-1 


where  Sn-1  is  a  function  that  depends  on  lower  order  terms 
q j<n.  It  is  easy  to  show  by  the  method  of  varia¬ 
tion  of  parameters  that 


f  =  C  f  (z) 
n  no 


4* 

-  f  1  n 


(n)  dn 


qn  =  -Cn  exp(a ( z )  )/E  +  ln(z) 


where 


=  f  g  (n 


)  dn 


then  formula  (34)  would  be 


Az . 


Az2 


=  ii_i  +  (fi  +  fi-i>  +  IF  (-fi'  +  "i-i*1 


Az3 


120 


I  <fi"  +  fi-lM)  * 


3.4 


RESULTS  FOR  THE  SADDLE  POINT 


The  series  expansion  method  outlined  in  the  pre¬ 
vious  sections  was  programmed.  A  calculation  was  made  for  50 
terms  of  this  series.  The  aim  was  to  apply  the  techniques 
expounded  by  Van  Dyke  (1984)  to  determine  the  nature  of  the 
true  solution  of  the  problem.  The  calculations  were  per¬ 
formed  once  in  single  precision  and  once  in  double  precision 
to  determine  the  effects  of  roundoff  error.  By  comparing  the 
results  of  these  two  calculations  it  was  found  that  roundoff 
errors  are  probably  confined  to  the  last  four  of  the  16 

decimal  digits  of  the  double  precision  results.  Thus,  the 

12 

accurate  numerical  integration  method  (one  part  in  10  )  and 
stable  roundoff  behavior  indicates  that  the  computed  series 
coefficients  may  be  accurate  to  10  or  12  decimal  places. 

It  is  not  effective  to  analyize  the  complete  ve¬ 
locity  field  (h,r,f,q).  Instead  the  wall  skin  friction  com¬ 
ponent  in  the  attachment  line  direction  ( x-direction)  is 
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considered  for  detailed  study.  The  wall  skin  friction  coef 


ficient  Cfx  is  given  by 


C _  —  T  /PU 

fx  wx  ° 


=  v/v  e  sinx  (3f/3z)  q/Uq 


(43) 


so  that 


U2  C,  / xTv  e  sin  =  -q(0)  (44) 

o  fx  v  x 

Thus,  the  series  expansion  for  q(0)  is  examined  in  detail. 
These  results  can  then  also  be  directly  compared  to  Davey's 
(1961)  results. 

The  series  for  q(Q)  which  will  be  denoted  simply  by 
t  for  simplicity  can  be  written  as 

N 

t  =  q (0)  =  ^  entn  (45) 

n=o 


where  N=50,  but  other,  lower,  values  of  N  will  also  be  con- 

The  complete  table  of  coefficients  t  for  N«50  is 

n 


sidered. 


given  in  Appendix  A.  These  coefficients  are  rounded  to  10 
figures. 


The  coefficients  t  increase  in  magnitude  and  have 
a  sign  pattern  that  indicates  complex  conjugate  singularities 
near  the  origin  of  the  complex  e-plane.  Any  complex  singu¬ 
larity  must  be  paired  with  its  complex  conjugate  because  the 

coefficients  t  are  all  real, 
n 


A  Domb-Sykes  plot  (see  Van  Dyke,  1984)  of  the 
coefficients  is  very  irregular  and  does  not  give  useful 
information.  A  similar  Domb-Sykes  type  plot  using  the  root 
test 


e*  =  limit  |  t  | 
n 


gives  an  estimate  of  the  radius  of  convergence  e*  of  this 
series.  In  a  Domb-Sykes  plot  the  quantity  1 1 n I  1^n  is  plot¬ 
ted  versus  values  of  1/n  and  extrapolated  linearly  to  (1/n)  ■ 
0  (i.e.,  n=«) .  The  value  so  obtained  is  an  estimate  of  et. 
It  was  found  in  this  way  that 


v$*a 


Thus,  the  radius  of  convergence  of  the  series  expansion  is 
about  0.55.  This  does  cover  a  range  of  values  of  e  for  which 
backflow  along  the  attachment  line  does  occur.  Recall  that 
Davey  (1961)  showed  that  backflow  at  the  saddle  point  (i.e., 
flow  away  from  the  saddle  point  along  the  attachment  line, 
whereas  the  outer  inviscid  flow  is  towards  the  saddle  point) 
occurs  for  e=0.4294.  However,  the  limit  e*=0.55  is  a  disap- 
pointly  small  fraction  of  the  range  e=[0,l]  for  which  the 
saddle  point  flow  exists. 


In  order  to  overcome  this  limitation  on  the  conver¬ 
gence  of  the  series  (45),  it  can  be  recast  in  terms  of  Pade' 
fractions  (Van  Dyke,  1984).  This  recasting  of  a  series 
effectively  serves  as  a  device  to  analytically  continue  t 
outside  the  region  of  convergence.  The  Pade'  fractions  cap¬ 
ture  the  singularities  that  limit  convergence  as  poles  in  the 
denominator.  Branch  point  singularities  are  approximated  by 
a  distribution  of  very  weak  poles  along  the  branch  cuts. 
These  properties  of  the  Pade'  fractions  emerge  automatically 
from  the  mere  calculation  of  the  coefficients. 


Thus,  the  series  (45)  is  written  as 


N 


n=o 


£  v1 


1=0 


M 


1  +  E 

j  =  l  ^ 


(47) 
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where  L+M=N.  By  reexpanding  the  Pade'  fraction  on  the  right 
of  equation  (47)  into  a  power  series  about  e*0  one  can  match 
the  coefficients  of  the  powers  of  e  and  obtain  a  linear 
system  of  equations  for  the  coefficient  a^  and  b^.  Once 
these  coefficients  are  in  hand  the  numerator  and  denominator 
are  solved  for  their  respective  zeroes  to  determine  the 
zeroes  and  poles  of  the  Pade*  fraction.  This  type  of  Pade' 
analysis  can  be  carried  out  for  all  combinations  of  0<L+M£N. 
In  most  cases  only  the  fractions  with  L=M  are  useful. 

In  the  present  study  the  following  Pade'  fractions 
were  analyized;  (L,M)  =  (K-1,K),  (K,K)  and  (K,K-1)  where 

11  <K£(2N+l)/2.  It  was  found  that  the  (K-1,K)  and  (K,K-1) 
fractions  behaved  erratically  and  seemed  to  converge  to  use¬ 
ful  results  only  for  values  of  e<0.7.  The  equal  order  frac¬ 
tion  (K,K)  produced  consistant  and  convergent  results  for 
values  of  e  fairly  close  to  the  limiting  value  e=l.  The 
(K,K)  Pade'  fraction  for  K=25,  corresponding  to  N=50,  is 
tabulated  in  Appendix  A. 

The  results  of  this  calculation  and  the  analysis  of 
the  Pade'  fraction  indicates  the  following  facts: 

(i)  The  function  t(e)  has  a  zero  (corresponding  to  the 
onset  of  reverse  flow)  of  e=0 . 429368226 .  This  number 
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seems  to  be  correct  based  on  various  orders  of  trunca¬ 


tion,  numerical  integration  step  size  and  placement  of 
the  outer  integration  limit  zm.  This  value  also  com¬ 
pares  favorably  to  Davey's  value  of  0.4294  and  the 
improved  value  of  Cooke  and  Robins  (1970)  of 
e=0. 42937. 

(ii)  The  singularities  blocking  the  convergence  of  the 
series  (45)  are  a  pair  of  branch  points  located  at 

e  =  0.427  +  i0.346 

where  i  =  n/~1.  Note  that  this  gives  a  radius  of 
convergence  e*  *  0.55  confirming  the  value  obtained  by 
the  root  test.  There  are  poles  all  along  two  arcs, 
one  from  each  of  the  two  branch  points,  extending 
towards  positive  values  of  Real(e).  The  residue  of 
each  of  these  poles  is  very  small.  hs  the  value  of  K 
is  increased  more  and  more  of  these  poles  appear,  thus 
leading  to  the  conclusion  that  the  above  singularities 
are  branch  points.  Unfortunately  the  irregular  be¬ 
havior  of  the  Domb-Sykes  plot  prevents  the  determina¬ 
tion  of  the  exponent  of  these  branch  points. 


(iii)  Since  the  branch  points  appear  for  Real(e)  >  0,  the 

Euler  transformation  that  maps  these  points  to 


infinity  is  not  useful.  Thus  no  information  can  be 
obtained  by  an  Euler  transformation.  However,  if  the 
series  values  of  e<0  were  desired  (i.e.,  the  nodal 
point  analysis)  then  the  Euler  transformation  would 
yield  a  useful  new  series. 

(iv)  The  Pade'  fraction  has  no  singularities  on  -1  £  Real 
(e)  <  1  hence  it  can  be  used  to  evaluate  t  for  real 

values  of  e.  However,  the  accuracy  of  these  values 
can  only  be  tested  experimentally  by  monitoring  the 
values  of  t  at  varying  orders  of  truncation  K.  Table 
II  gives  the  calculated  values  of  t  as  a  function  of  e 
for  K=10 ( N=20 ) ,  K=15 { N=30 ) ,  K=25(N=50)  and  Davey's 

values.  It  can  be  seen  from  Table  II  that  convergence 
of  the  values  of  t  is  very  rapid  in  K  for  values  of 
e£0.6.  However,  for  values  of  e>0.6  the  convergence 
rate  with  respect  to  K  is  fairly  slow  especially  for 
e>0.8.  For  e=0.9,  K=25  (N=50  terms)  produces  only  two 
significant  figures  of  accuracy.  It  is  important  to 
note  that  for  K»10(N=20)  the  series  expansion  method 
will  give  reliable  values  of  t  only  up  to  e=0.7.  For 
the  general  case  of  the  attachment  line  problem  the 
series  expansion  method  is  a  double  series  in  x  and  e 
so  that  calculating  a  large  number  of  terms,  such  as 
N=50  may  not  be  feasible.  Nevertheless,  the  method 
can  give  accurate  results  with  K*10  at  least  for  e<0.7 

3-18 


A/./ 


vlv**-s 


-Vv' 


«■  / 


which  is  substantially  within  the  backflow  region 
e>0 . 4294 . 

It  should  be  mentioned  that  we  have  also  solved 
this  problem  by  direct  numerical  integration  of  the  saddle 
point  equations.  This  calculation  confirms  Davey's  results 
for  all  his  values  of  e  to  the  four  significant  figures  that 
he  gives. 

In  summary,  it  should  be  noted  that  the  series 
expansion  method  does  provide  a  useful  way  of  evaluating  the 
saddle  point  flow  for  values  of  e  at  least  up  to  0.7  with  a 
modest  number  of  terms.  The  series  convergence  is  blocked  by 
branch  points  at  a  radius  of  about  0.55  and  these  cannot 
be  removed  easily. 

In  Figure  3  the  velocity  profile  f(z)  at  the  saddle 
point  for  e=0.7  is  plotted.  The  profile  shows  the  reversed 
flow  region  at  the  wall  and  shows  that  separation  occurs  on 
the  attachment  line  ahead  of  the  saddle  point.  The  question 
that  one  faces  now  is  whether  the  saddle  point  solutions  for 
e>0.4294  have  any  relevance  when  separation  occurs  upstream. 
This  question  is  to  be  answered  by  analyizing  the  entire 
attachment  line  boundary  layer. 


Section  4 


THE  PLOW  ALONG  THE  ATTACHMENT  LINE-DIRECT 
NUMERICAL  SOLUTION 

4.1  DESCRIPTION  OF  THE  FLOW  AND  NUMERICAL  METHOD 

The  aim  now  is  to  solve  the  complete  attachment 
line  problem  formulated  by  equations  {12  -  18).  The  primary 
question  which  we  wish  to  answer  is  the  following:  Can  the 

boundary  layer  proceed  through  the  separation  point,  which 
occurs  upstream  of  the  saddle  point,  and  terminates  at  the 
saddle  point  with  Davey's  separated  solution?  Implicit  in 
this  question  is  if  this  separated  boundary  layer  profoundly 
affects  the  external  inviscid  flow  (i.e.,  5  a-*00  at  separa¬ 

tion)  and  viscous/inviscid  interactions  are  required.  If  the 
outer  potential  flow  does  require  modifications  in  order  for 
the  boundary  layer  to  proceed  past  separation,  then  Davey's 
solutions  for  e>0.4294  are  irrelevant.  Davey  (1961)  and 
Cooke  and  Robins  (1970)  assumed  this  to  be  so.  In  parti¬ 
cular,  Cooke  and  Robins  (1970)  examined  this  problem,  but  for 
an  artificial  inviscid  flow,  for  e<0.4294  and  simply  re¬ 
iterated  Davey's  assumption  that  for  e>0.4294  the  boundary 
layer  cannot  proceed  past  the  separation  point  on  the  at¬ 
tachment  line. 


It  is  not  at  all  obvious  that  Davey's  assertion  is 

correct.  Separation  can  occur  with  the  displacement 

thickness  S*  remaining  bounded.  Under  such  a  circumstance 

-1/2 

the  separated  flow  does  not,  to  order  R  '  where  R  is 
Reynolds  number,  disturb  the  inviscid  flow  (Klemp  and 
Acrivos,  1976).  It  is  a  common  occurance  in  three  dimen¬ 
sional  symmetry  plane  boundary  layers  that  the  displacement 
thickness  not  only  remains  bounded  but  in  fact  decreases  as 
separation  is  approached.  An  example  of  this  is  given  in 
Part  II  of  this  report  concerning  the  leading  edge  boundary 
layer  on  a  disk  in  edgewise  flow.  The  feature  that  causes 
the  boundary  layer  thickness  to  remain  bounded  even  at 
separation  is  the  very  strong  lateral  outflow  to  either  side 
of  the  symmetry  plane.  These  flow  features  are  present  in 
the  attachment  line  boundary  layer  considered  here. 

The  equations  (12  -  18)  are  solved  numerically  by 
starting  at  the  nodal  point  x=0  and  marching  the  solution  in 
increasing  values  of  x  until  the  saddle  point  x®77  is  reached. 
The  numerical  method  used  is  the  standard  box  method  and 
needs  no  further  elaboration.  Reference  can  be  made  to 
Cebeci  and  Bradshaw  (1977)  for  the  details  of  applying  the 
box  method  to  general  boundary  layer  equations. 

It  should  be  noted  that  for  values  of  e>0.4294 
separation  of  the  f  component  of  velocity  will  occur  ahead  of 


the  saddle  point  x=tt  and  thus  the  numerical  method  must  march 
against  the  stream,  f<0,  in  part  of  the  boundary  layer.  This 
problem  can  be  circumvented  by  noting  that  the  cases  for 
which  f<0  have  f  <<  1  so  that  the  terms 


in  equations  (14)  and  (15)  can  be  set  to  zero  when  f£0.  This 
avoids  the  unstable  phenomena  of  effectively  negative  dif¬ 
fusion  (i.e.,  marching  upstream  instead  of  downstream).  Such 
a  procedure  is  now  called  the  FLARE  approximation.  Its 
accuracy  and  general  validity  is  not  completely  quantified  so 
that  results  obtained  using  this  approximation  must  be  ac¬ 
cepted  with  caution.  Also  the  FLARE  approximation  may  not 
work  for  flow  velocities  f<0  whose  magnitude  exceeds  a 
certain  level.  For  these  reasons  it  is  necessary  to  calcu¬ 
late  the  flow  with  e>0.4294  with  other  methods  that  do  not 
make  use  of  such  a  crude  approximation  and  that  are  more 
robust  for  larger  values  of  e.  This  motivated  our  work  to 
solve  equations  (12  -  18)  by  high  order  perturbation  series 
in  e. 

4.2  SOME  RESULTS  OF  THE  NUMERICAL  INTEGRATION 

The  direct  numerical  integration  method  was  applied 
to  two  cases  of  nodal  to  saddle  attachment  line  flow.  The 
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first  case  was  to  the  flow  described  by  the  inviscid  flow  of 


equation  (1).  This  calculation  resulted  in  successful  inte¬ 
grations  from  the  nodal  to  the  saddle  point  only  for  values 
of  e£0.55.  This  range  of  values  of  e  does  encompass  the 
separated  flow,  but  it  was  felt  that  the  complete  boundary 
layer  can  be  calculated  for  larger  values  of  e.  The  value  of 
0.55  at  which  the  FLARE  approximation  ceases  to  work  has 
nothing  to  do  with  the  value  of  e*=0.55  for  the  convergence 
of  the  saddle  point  series. 

The  second  case  considered  was  the  inviscid  flow 
used  by  Cooke  and  Robins  (1970).  Their  inviscid  flow  also 
corresponded  to  an  attachment  line  flow  between  nodal  and 
saddle  points,  but  it  did  not  correspond  to  any  particular 
geometric  shape  and  thus  was  somewhat  artificial.  However  it 
was  useful  to  compare  our  calculations  with  theirs,  so  their 
case  was  calculated.  Since  the  results  and  conclusions  are 
qualitatively  the  same  for  either  external  inviscid  flow,  the 
case  of  Cooke  and  Robins  will  suffice  to  illustrate  points. 

The  Cooke  and  Robins  case  was  prepared  as  a  short 
paper  and  submitted  to  the  Journal  of  Applied  Mechanics. 
This  paper  is  reproduced  in  Appendix  B.  Thus,  one  can  turn 
to  this  appendix  for  the  details  of  the  results. 


It  is  sufficient  to  state  here  that  the  boundary 
layer  on  the  attachment  line  can  indeed  be  computed  through 
separation  all  the  way  to  the  saddle  point  for  values  of  e 
larger  than  0.4294.  (Notes  Since  the  external  inviscid  flow 
has  different  forms  in  the  present  and  Cooke  and  Robins 
cases,  the  values  of  e  do  not  correspond  one-to-one  on  the 
attachment  line.  However,  at  the  saddle  point  the  values  of 
e  do  correspond.  Thus  e>0.4294  always  implies  separation 
ahead  of  the  saddle  point,  but  at  different  locations  for 
different  external  inviscid  flows.) 

As  mentioned  in  Appendix  B,  it  seems  reasonable  to 
expect  that  the  boundary  layer  flow  terminating  with  Davey's 
saddle  point  solution  exists  for  even  larger  values  of  e. 
For  this  reason  the  e-expansion  techniques  was  used  in  the 


next  section. 


Section  5 


PERTURBATION  SERIES  SOLUTION  OF  THE  ATTACHMENT  LINE  FLOW 

5.1  BASIC  FORMULATION  AND  EXPANSION 

In  order  to  both  verify  and  extend  the  solutions 
obtained  by  the  direct  numerical  integration  of  Section  4, 
the  attachment  line  equations  are  now  solved  by  a  series 
expansion  in  both  x  and  e.  The  starting  point  is  equations 
(12  -  18). 


First  it  is  convenient  to  make  the  change  in  inde¬ 
pendent  variable  from  x  to  s  by  the  transformation 

s  =  -  cos  x  .  (48) 


Equations  (12  -  16)  then  become 


e(1-s  >  f  9¥  +  al  -  gr  +  n 


± 


149C) 


ti  +  <  = 


o 


( 49d) 


e  ^ ( 1-s2) 


es 


( 49e ) 


The  boundary  conditions  remain  equations  (17  -  18).  Next  the 
vector  F  =  (h,g,r,f,q)+  is  expanded  as 

*’t  t eV" f-  •  <5<” 

n=o  m=o 

By  substituting  expansion  (50)  into  equations  (49)  and  col¬ 
lecting  the  appropriate  terms  of  equal  powers  in  ensm  one 
obtains  the  following  perturbation  equations: 


—iiili  -  (g  r  +  r  g  -2hh  '  =  P 
9z  o  nm  o  nm  o  nm  rcnm 


+  q  =0 
9z  nm 


9z  ^o^nm  ^enm 


where  the  prime  on  the  parenthesis  in  equation  (5L 
indicates  that  the  term  is  multiplied  by  1/2  for  n=m=0.  T1 
functions  Pbnm»  Pcnm  and  Penm  ace  define<J  as  follows: 


p,  =  m  f  ,  i  -  (m+1)  f  . 
bnm  n-l,m-l  n-l,m+l 


n-1  l 


p  =  66  + 

^cnm  no  mo 


Z)  £  (g  krn-*.m-k 


1=0  k=o 
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-  (m-k)  fn_t-l,n,-k- 

where 

{1  for  m  <  n-1 
0  for  m  >_  n-1 

The  details  of  the  derivation  of  these  equations 
are  straightforward,  but  somewhat  tedious.  The  main  points 
to  note  are  the  following: 

(i)  Equations  51  and  52  hold  for  all  values  of  n>0  and  for 
0<m<n  for  each  n.  Furthermore  g^^  =  g^  etc.  If  the 

upper  limit  in  one  of  the  sums  in  Eqs.  (52)  is  less 
than  the  lower  limit  then  the  sum  is  null.  <5^  is  the 


kronecker  delta  function. 


(ii)  For  n+m  =  odd  F _ =0.  This  greatly  reduces  the  amount 

nm 

of  storage  required  for  large  values  of  N.  The  amount 
of  computer  storage  available  is  the  main  limiting 
factor  on  the  number  N  of  terms  that  can  be  computed. 

(iii)  The  differential  operators  on  the  left  side  of  the 

equations  (51)  are  identical  to  the  differential 
operators  of  the  saddle  point  expansion  of  Section  3. 
Only  the  right  hand  sides  of  equations  (51)  are  dif¬ 
ferent.  Thus  the  exact  same  numerical  techniques  and 
also  about  75%  of  the  computer  program  for  the  saddle 
point  expansion  can  now  be  carried  over  to  the  at¬ 

tachment  line  equation. 

->■  ■> 

(iv)  The  boundary  conditions  on  FQo  =  Fq  are  the  same  as 

Eqs .  (28)  and  those  for  F„.  n>0  are  the  same  as  Eqs. 

nm 

(30)  . 


The  numerical  procedure  for  calculating  the  coef¬ 


ficients  F  need  no  further  discussion  other  than  to  note 
nm 


that  there  are  many  more  terms.  In  fact  the  number 


of 


terms  F _  is  given  by  the  formula 

nm 


r  T(n+d2i 
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rN+11 

4  4  J 
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for  N  even 


for  N  odd 


Thus,  for  N»20,  T=121  terms  are  required.  Since  Pnm  is  a  5- 
dimensional  vector,  and  if  100  points  in  z  are  retained 
across  the  interval  [0,zm]  then  a  total  of  60500  double 
precision  numbers  must  be  stored.  But  in  order  to  keep  the 
program  complexity  and  computing  times  reasonable  certain  of 
the  components  of  F*  and  F''  „  must  also  be  stored.  This 
basically  doubles  the  above  number  of  double  precision 
numbers  to  120000.  Since  one  double  precision  number 
requires  8  bytes  of  storage  a  total  of  about  1  million  bytes 
of  storage  is  required.  This  amount  of  storage  is  about  the 
limit  available  on  our  machine.  Thus,  the  present  calcula¬ 
tion  is  limited  to  N=20.  However,  it  is  a  simple  matter  to 
transport  our  program  to  a  larger  machine  and  obtain  results 
for  larger  values  of  N  if  such  a  machine  becomes  available. 
The  computer  time  required  for  N-20  is  only  15  minutes  on  a 
DEC-10. 


5.2  RESULTS  OP  THE  PERTURBATION  EXPANSION  OP  THE 

ATTACHMENT  LINE  PLOW 


A  total  of  N»20  terms  (actually  21  terms  including 
n®0)  were  computed.  Again  we  will  illustrate  our  results  by 


simply  considering  the  value  of  (3f/3z)_  n  ■  -q(o)  which  is 

z*u 

proportional  to  the  skin  friction  in  the  attachment  line 
direction.  Thus  the  series 


q  (0) 


n  n 


-E  £ 


n  m 
e  s 


(53) 


is  described.  It  should  be  noted  that  the  attachment  line 
direction  skin  friction  coefficient  Cfx  =  eyfv  sinxq(0 )/UQ2. 

The  series  (53)  was  first  summed  for  each  value  of 

s  as 


n 

V01  =  2  qnm(0) 

m=o 


to  give  the  series 


(54) 


N 

q  ( 0)  =  ^  e11  qn(0)  . 

n=o 


(55) 


This  later  series  is  then  recast  as  a  Pade'  fraction  and 


evaluated  for  various  values  of  e 


It  was  found  that  for  e<0.55  the  series  was  conver¬ 
gent  for  all  values  of  s  and  gave  essentially  identical 
results  to  the  direct  numerical  method  of  Section  4,  even  in 
the  separated  region.  For  e<0.55,  the  Pade'  analysis  was 
unnecessary  since  the  series  converges  for  all  values  of  s. 
These#  relatively  uninteresting  results  will  thus  not  be 
shown. 

For  values  of  e>0.55  it  was  already  indicated  that 
the  direct  numerical  integration  method  could  not  proceed 
beyond  the  point  of  separation.  However/  the  series  expan¬ 
sion  method  with  N=20  and  with  the  help  of  the  (K,K)  Pade' 
fraction  with  K=10  did  seem  to  produce  accurate  values  of 
q( 0 )  for  all  values  of  sin  [-1,1]  i.e.,  x  in  [0,tt]  and  for 
e=0.6  and  possibly  also  e=0.65.  Figure  4  shows  the  plot  of 
q(0)  versus  x  for  e=0.6,  0.65  and  0.7  obtained  by  the  series 
method  and  the  numerical  integration.  In  all  three  cases  the 
direct  numerical  integration  yields  values  of  q(0)  which 
nearly  coincide  with  the  values  of  q(0)  obtained  by  the 
series  up  to  the  point  of  separation.  Even  for  e*0.7, 
separation  is  predicted  at  x=0.71  by  the  direct  numerical 
method  and  x=0.714  by  the  series  method.  However,  beyond 
separation  the  predicted  values  of  q(0)  for  e»0.7  behave  in 
an  irregular  manner  and  thus  leave  some  doubt  as  to  the 
convergence  of  the  series  method.  It  was  already  seen  in  the 
saddle  point  expansion  that  N*20  is  inadequate  to  get  more 


than  possibly  two  significant  figures  for  q{0).  There  seems 
to  be  some  indication  in  Figure  3  that  the  convergence  of  the 
Pade'  approximation  is  poorer  just  beyond  separation  than  at 
the  saddle  point.  If  this  is  so,  then  there  is  a  possibility 
that  the  flow  does  become  singular  at  separation  for  values 
of  e<1.0  and  cannot  be  continued,  without  a  modification  to 
the  external  inviscid  flow,  all  the  way  to  the  saddle  point. 
In  this  case  Davey's  saddle  point  solutions  are  irrelevant. 
However,  it  does  seem  clear  that  the  flow  on  the  attachment 
line  does  pass  smoothly  through  separation  and  terminates  at 
the  saddle  point  with  Davey's  profile  for  values  of  e  at 
least  up  to  0.65.  Thus  there  is  a  substantial  range  of  flows 
exhibiting  the  separated  horseshoe  vortex  that  can  be  calcu¬ 
lated  by  pure  boundary  layer  theory  with  no  viscous/inviscid 
interactions. 

Figure  5  shows  a  sequence  of  attachment  line  direc¬ 
tion  velocity  profiles  in  the  neighborhood  of  the  saddle 
point  for  e=0.55.  These  profiles  were  calculated  using  the 
direct  numerical  integration.  The  figure  clearly  shows  the 
small  separated  region.  Since  the  attachment  line  is  in  the 
spanwise  direction  of  the  cylinder  and  the  main  flow  direc¬ 
tion  is  around  the  circumference,  the  attachment  line  flow  is 
actually  the  start  of  the  cross-flow  of  the  main  three- 
dimensional  boundary  layer  of  the  cylinder.  Also  since  the 
vorticity  generated  by  the  attachment  line  component  of  the 


velocity  is  perpendicular  to  the  attachment  line.  Figure  5 
indicates  the  existence  of  a  streamwise  (circumferential) 
concentrated  vortex  that  will  surround  the  cylinder  at  the 
narrow  positions. 


5-10 


CONCLUDING  REMARKS 


The  calculation  described  in  the  previous  sections 
have  shown  the  interesting  fact  that  along  the  line  of  at¬ 
tachment  on  the  plane  of  symmetry  of  three  dimensional  flows 
it  is  possible  to  have  a  separated  boundary  layer  that  is 
determined  by  the  outer  inviscid  stream  corresponding  to  the 
body  alone.  This  is  an  important  finding  because  it  indi¬ 
cates  that  in  some  cases  the  horseshoe  vortex  formed  by  the 
boundary  layer  rollup  ahead  of  an  obstruction  can  still  be 
calculated  by  boundary  layer  theory  alone.  A  Navier-Stokes, 
or  even  a  PNS,  calculation  is  not  required  to  predict  this 
horseshoe  vortex. 

It  was  also  learned  that  the  FLARE  method  for 
computing  in  the  separated  region  is  inadequate  if  the 
separation  is  too  large,  even  though  the  boundary  layer 
theory  still  holds.  The  perturbation  method  is  not  a  suit¬ 
able  computational  tool  to  deal  with  practical  geometries 
such  as  a  destroyer  bow.  Such  shapes  do  not  lend  themselves 
to  simple  one  dimensional  parameter izations.  However  a 
practical  method  for  handling  separated  flow  or  the  general 
geometries  might  be  the  method  of  Klemp  and  Acrivas  (1976). 
This  certainly  should  be  investigated.  It  is  important  to 


develop  such  calculation  methods  so  that  general  destroyer 
bow  domes  can  be  treated.  This  is  a  necessary  step  if  a 
rational  boundary  layer  procedure  is  required  for  destroyer 
hulls. 

These  same  comments  hold  for  the  case  of  the  root 
boundary  layer  of  the  blade/hub  intersection.  In  this  prob¬ 
lem  it  can  be  seen  that  a  certain  amount  of  filleting  can 
reduce  the  intersection  horseshoe  vortex  to  such  a  size  that 
it  does  not  disturb  in  a  gross  way  the  body  generated  po¬ 
tential  flow.  This  is  all  that  is  really  needed  for  good 
fillet  design. 

On  full  scale  ships  the  region  of  laminar  flow  is 
likely  to  be  a  tiny  part  of  the  hull  surrounding  the  nodal 
point  of  attachment.  Surely  the  boundary  layer  will  transi¬ 
tion  far  ahead  of  the  saddle  point.  There  is  no  reason  to 
believe  that  the  attachment  line  boundary  layer  could  not  as 
easily  be  calculated  for  turbulent  flow  with  a  suitable 
turbulence  model.  However,  since  the  velocity  decreases  as 
the  saddle  point  of  attachment  is  approached  and  since  the 
divergence  of  the  flow  is  very  strong  there  it  may  happen 
that  the  flow  will  retransition  back  to  laminar  flow  before 
the  saddle  point  is  reached.  This  is  an  interesting  question 
that  needs  further  investigation  and  would  have  serious 
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Figure  2.  The  flow  attachment  at  the  bow  of  a  destroyer  hull. 

Point  N  is  a  nodal  point  and  S  is  a  saddle  point  of 
attachment . 


The  velocity  profiles  u(z)  for  the  flow  along  the 
attachment  line  for  e=0.55.  The  saddle  point  is 
at  x=ir. 


Appendix  A 


TABLES  OF  NUMERICAL  RESULTS 


SERIES (50) (X) 


PADE(25, 25) (X)  = 


- . 5704652525D+00*X**  0 
. 1053234471D+01*X**  1 
. 5398605816D+00*X**  2 
. 3632211708D+00*X**  3 
. 1048701481D+00*X**  4 
- . 3171103056D+00*X**  5 
>.8735167146D+00*X**  6 
>.1353168116D+01*X**  7 
1288116911D+01*X**  8 
. 4580482502D_02*X**  9 
. 3094025117D+01*X**10 
. 7660909416D+01*X**11 
. 1132759345D+02*X**12 
. 868626931 2D+01*X**1 3 
- . 8003479724D+01*X**14 
- . 4416612966D+02*X**15 
- . 9239648057D+02*X**16 
1175639805P+03*X**17 
- . 4721269869D+02*X**18 
. 2107858533D+03*X**19 
. 6922902362D+03*X*#20 
. 1223617160D+04*X**21 
. 1232028289D+04*X**22 
- . 3106514888D+03*X**23 
- . 4472868238D+04*X**24 
- . 1106525063D+05*X**25 
- . 1 63624229 1D+05*X**26 
- . 108352067 1D+05*X**27 
. 2010132069D+05*X**28 
. 8712059342D+05*X**29 
. 1747726846D+06*X**30 
. 2082729B18D+06*X**31 
. 3098506788D+05*X**32 
5547394114D+06*X**33 
- . 1606875604D+07*X**34 
- . 2666760147D+07*X**35 
- . 2312212941D+07*X**36 
. 1949596261D+07*X**37 
. 1251828419D+08*X**38 
. 2820746620D+08*X**39 
. 3825127696D+08*X**40 
. 1674802350D+08*X**41 
- . 7382586266D+08*X**42 
- . 2554261 145D+09*X**4 3 
4704391168D+09*X**44 
4906990721D+09*X**45 
. 1316733770D+09*X**46 
. 1912208305D+10*X**47 
. 485746020 4D+10*X**48 
. 7351166020D+10*X**49 
. 4901180938D+10*X**50 


- . 6704652525D+00*X**  0 
. 7059840489D+01*X**  1 
3739315631D+02*X**  2 
. 1017958962D+03*X*»  3 
-.1083245618D+03*X**  4 
1660314591D+03*X**  5 
. 6932322542D+03*X**  6 
- . 7877205048D+03#X**  7 
. 1804664342D+03*X*»  8 
7824261432D+03*X*#  9 
. 4189294352D+04*X**10 
- . 6161979200D+04*X**11 
. 1459220753D+04*X**12 
. 2717835772D+04*X**13 
. 3322742735D+04*X**14 
. 1127274426D+04*X**15 
5539762684D+05*X**18 
. 1481042726D+06*X**17 
- . 2032932536D+06*X**18 
. 1715036965D+06*X**19 
- . 9106075210D+05*X**20 
. 2848786425D+05*X**21 
- . 4053012796D+04*X**22 
- . 1209159580D+03*X**23 
. 8174016797D+02*X**24 
38768B7074D+01*X**25 


.1000000000D+01*X**  0 
1052931093D+02*X**  1 
. 4705489693D+02*X**  2 
- . 1008951820D+03*X**  3 
. 4161824668D+02*X**  4 
. 2998708736D+03*X**  5 
- . 6734451159D+03*X**  6 
,4168034093D+03*X**  7 
. 2120067392D+02#X**  8 
. 1475005815D+04*X**  9 
- . 4550740078D+04*X**10 
. 3886743493D+04*X**ll 
. 1381576460D+04*X**12 
1000111811D+04*X**13 
- . 5412715045D+04*X**14 
- . 1107880567D+05*X**15 
. 7207196025D+05*X**16 
1408489129D+06*X**17 
. 1545824681D+06*X**18 
1054640933D+06*X**19 
. 4410557307D+05*X**20 
- . 9943530505D+04*X**21 
. 6315803747D+03*X**22 
. 1426783601D+03*X**23 
- . 1B56963426D+02*X**24 
. 2515278745D+00*X**25 


Table  I 


The  saddle  point  expansion  coefficients  for  N=50  series  and  the 
(25,25)  Pade’  expansion.  Here  X  stands  for  the  expansion  parame 
ter  e  in  the  text. 


K=10 

K=  15 

K=25 

Davey 

57046525 

.57046525 

.57946525 

0.5705 

45937368 

.45937368 

.45937368 

0.4594 

33532784 

.33532784 

.33532784 

0.3353 

19700104 

.19700104 

.19700104 

0.1970 

04596812 

.04596814 

.04596814 

0.0460 

11150569 

-.11150066 

-.11150066 

-0.1115 

26689931 

-.26659399 

-.26659376 

-0.2666 

41834578 

-.41301722 

-.41297965 

-0.4130 

58709785 

-.55046973 

-.54931101 

-0.5488 

82489666 

-.69458402 

-.68365346 

-0.6761 

12049790 

-.89837162 

-.85222725 

-0.8112 

Table  II 


The  values  of  -q(0)  at  the  saddle  point  of  attachment  for 
various  truncations  N  of  the  series  expansion  and  the  (K,K) 
Pade'  approximation,  K=N/2,  compared  to  Davey's  numerical 
results. 


THE  BOUNDARY  LAYER  ON  AN  OBLATE 
SPHERIOD  AT  AN  ANGLE  OF  ATTACK 


Section  1 


INTRODUCTION 

The  aim  of  this  research  project  was  to  investigate 
leading  edge  boundary  layers  on  wings  and  propeller  blades. 
The  reason  for  this  investigation  is  that  existing  three 
dimensional  boundary  layer  methods  (see  Groves,  1981,  Groves 
and  Chang,  1985,  Cebeci,  1985)  do  not  account  properly  for 
such  leading  edge  boundary  layers.  For  example.  Groves 
(1981)  simply  neglects  the  leading  edge  and  begins  the  calcu¬ 
lation  on  the  face  of  the  blade  with  flat  plate  turbulent 
boundary  layer  data.  Groves  and  Chang  (1985)  assume  that  the 
leading  edge  is  wedge  shaped  and  start  their  calculation  with 
basically  the  Falkner-Skin  similarity  profiles.  However,  for 
subsonic  flows  wing  and  blade  leading  edges  are  rounded,  not 
sharp.  Cebeci  (1985)  mentions,  in  his  review  of  wing 
boundary  layers,  only  the  case  where  the  leading  edge  is  a 
symmetry  plane.  In  this  case  the  leading  edge  is  an  at¬ 
tachment  line  along  which  the  crossflow  pressure  gradient 
vanishes.  This  reduces  the  boundary  layer  equations  to 
quasi-three  dimensional  form  that  can  be  solved  separately 
(as  initial  conditions)  from  the  remaining  three  dimensional 
equations  over  the  rest  of  the  surface.  Cebeci  (1985)  al¬ 
ludes  to  the  assumption  that  for  wings  this  type  of  leading 
edge  condition  is  approximately  true  even  for  small  angles  of 


attack.  This  assumption  is  probably  adequate  for  fairly  high 
aspect  ratio  wings  with  straight  leading  edges.  However,  for 
low  aspect  ratio  wings  at  moderate  angles  of  attack,  for 
marine  propellers  and  even  for  the  vicinity  of  the  wing  tip 
of  high  aspect  ratio  wings  the  attachament  line  curves  sub¬ 
stantially  off  the  leading  edge  and  is  not  a  symmetry  plane. 
Thus,  Cebeci's  assumption  would  not  be  valid  in  these  situa¬ 
tions. 

The  aim  was  in  fact  to  examine  the  leading  edge 
flow  in  the  vicinity  of  a  wing  or  blade  tip.  This  leading 
edge  flow  is  what  initiates  the  rollup  of  the  tip  vortex.  It 
seems  appropriate  then  to  study  the  leading  edge  boundary 
layer  on  a  flat  oblate  spheroid  in  edgewise  flow  as  depicted 
in  Figure  1.  (Figure  1  shows  a  rather  fat  oblate  spheroid 
only  for  the  purpose  of  more  clearly  depicting  some  of  the 
details.)  This  configuration  was  chosen  so  that  the  boundary 
layer  emanating  from  the  front  stagnation  point  would  proceed 
directly  to  the  rounded  tip  area  without  having  to  pass  over 
a  large  span  of  basically  two-dimensional  flow  that  is  rela¬ 
tively  benign  as  concerns  the  phenomena  under  consideration 
here . 

A  second  reason  for  choosing  the  oblate  spheroid  is 
that  its  geometry  and  potential  flow  is  analytically 


described.  This  allows  us  to  concentrate  on  the  boundary 
layer  theory  without  having  to  deal  with  complicated  numeri¬ 
cal  geometry  and  potential  flow  codes.  Figure  1  shows  an 
example  of  an  oblate  spheroid  in  relation  to  the  oncoming 
stream  whose  magnitude  is  and  angle  of  attack  is  a^.  Also 
a  surface,  nonorthogonal,  coordinate  system  used  for  boundary 
layer  calculations  is  shown.  In  the  course  of  this  investi¬ 
gation,  several  different  surface  coordinate  systems  had  to 
be  used.  All  these  are  easily  described  in  terms  of  the 
fundamental  elliptic  coordinate  describing  the  spheroid  and 
its  potential  flow. 


The  following  plan  describes  the  contents  of  this 
report.  section  2  contains  the  basic  spheroid  geometry, 
potential  flow,  various  surface  coordinate  systems,  the 
formulation  of  the  three-dimensional  boundary  layer  equations 
in  nonorthogonal  surface  coordinates  and  the  numerical  pro¬ 
cedures  for  solving  the  boundary  layer  equations.  Section  3 
describes  how  boundary  layer  calculation  methods  must  start. 
Sections  4  and  5  describe  and  discuss  the  results  obtained. 
Section  6  contains  some  concluding  remarks  and  prognosis  for 
continued  effort  on  the  problem. 


Section  2 
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THE  GEOMETRY  AND  BASIC  EQUATIONS  FOR  FLOW 
OVER  AN  OBLATE  SPINNING  SPHEROID 


2.1  GEOMETRY  AND  POTENTIAL  FLOW 


The  geometry  of  a  spheroid/  or  any  other  body#  is 
most  easily  described  in  terms  of  parametric  coordinates 
referred  back  to  a  global  cartesion  reference  frame.  Figure 
1  shows  an  oblate  spheroid  with  the  global  cartesian  (X/y,z) 
reference  frame  placed  at  its  centroid.  Then  the  surface  of 
this  spheroid,  whose  half-thickness  along  the  x-axis  is  r  and 
whose  radius  in  the  y-z  plane  is  1,  can  be  described  in 
vector  form  as 


R(y ,w) 
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(sinco]  +  cosojk) 


(1) 
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where  (i,j,k)  are  the  unit  vectors  corresponding  to  (x,y,z), 
(m,w)  are  the  or'hogonal  elliptic  coordinates  on  the  surface 
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of  the  spheroid  and  R(u,oj)  is  the  position  vector  from  the 


origin  of  the  (x,y,z)  frame  to  a  point  (u,w)  on  the  surface. 
The  length  scale  of  this  spheroid  is  the  radius  in  tne  y-z 
plane  which  has  been  normalized  to  1. 


The  potential  flow  whose  onset  velocity  is  in 
the  x-z  plane  can  be  found  in  Durand  Vol.  1  (1963).  Only  the 
surfaces  velocities  are  needed  in  this  work  so  that  only  the 
reduced  formulas  are  given  here.  Thus  the  surface  velocity 
is 
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where 


U  =  (y'/l-y2  -  a'4-M2)  /42+xV 


2  2 

U  =  Y  ,  Y  =  sin  a)  ,  X  =l=x 


(3) 


=  [xi  -  u(sinto3  +  cosujk)  /vl-jj  ' ]/h^ 


(4) 


e  =  coswj  -  sincuk 


h  =  ’/r2+X2p/'/i 
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h  =  vG.-y2/cosa) 
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and  a  is  proportional  to  the  angle  of  attack, 
the  angle  of  attack  defined  as 


In  terms  of 


(6) 
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the  value  of  a  is 


a  =  2a*[X<2-x2)  -  xcot-1  (|)] /[Xx-cot-1  (J)] 


(7) 


The  choice  of  the  proper  surface  coordinate  system 
in  which  to  perform  boundary  layer  calculations  is  a  delicate 
issue.  The  choice  depends  not  only  on  the  global  but  also  on 
the  local  topology  of  the  surface.  For  example,  in  the 
vicinity  of  the  stagnation  point  of  attachment  it  is  easiest 
to  deal  with  a  local  polar  coordinate  system  in  most  cases. 
This  is  particularly  true  for  bodies  whose  topology  near  the 
stagnation  point  is  nearly  spherical,  i.e.,  the  curvature  in 
every  direction  is  about  the  same  order  of  magnitude.  How¬ 
ever,  as  was  discovered  in  this  work,  for  a  flattened  body 
such  as  a  blade  or  wing,  or  disk,  such  a  system  is  not 
suitable.  Instead  a  locally  rectangular  system  needs  to  be 
used. 


In  this  work  the  basic  elliptic  system  was  not  used 
because  of  the  pole  at  the  top  and  bottom  (x=+x)  of  the  disk. 
Instead  two  other  surface  systems  were  used.  One  system  put 
the  poles  at  the  nose  and  tail  (z=+l)  at  the  edge  of  the  disk 


and  the  other  system  put  the  poles  at  the  tips  (y=+l).  We 
shall  refer  to  these  coordinate  systems  as  and  C2  respec¬ 
tively.  In  order  to  use  the  system  when  the  flow  is  at  an 
angle  of  attack,  which  moves  the  front  stagnation  point  off 
the  pole,  a  different  transformation,  given  by  Cebeci, 
Khattab  and  Stewartson  (1981)  was  used.  This  combined  system 

will  be  referred  to  as  C,  . 

lc 

It  is  best  to  leave  the  description  of  these  co¬ 
ordinate  systems  to  sections  where  they  are  applied.  At 
present  the  general  three  dimensional  laminar  boundary  layer 
equations  will  be  given  for  general  nonorthogonal  coordinate 
systems.  In  fact,  the  entire  development  and  computer  coding 
in  this  research  effort  was  carried  out  for  the  general  case, 
so  that  other  coordinate  systems,  body  geometries  and  po¬ 
tential  flows  could  be  easily  incorporated. 

2.2  THE  BOUNDARY  LAYER  EQUATIONS 

The  general  three-dimensional  boundary  layer  equa¬ 
tions  in  nonorthogonal  coordinates  can  be  found  in  several 
sources,  for  example,  Squire  (1955).  Let  (r,t,n)  be  such  a 
general  nonorthogonal  surface  coordinate  system.  The  (r,t) 
coordinates  lie  in  the  surface  and  n  can  be  assumed  to  be  an 
arclength  coordinate  tangent  to  the  body  normal.  Then  the 


metric  equation  for  an  element  ds  of  arclength  in  this  co¬ 
ordinate  system  is 

ds2  =  hr2  dr2  +  hfc2  dt2  +  2gdrdt  +  dn2  .  (8) 

The  position  R  of  any  point  on  the  body  surface  referred  to 
the  origin  of  the  global  cartesian  reference  frame  is  a 
function  of  (r,t)  and  can  be  written  as 

R  =  R(r ,t)  =  x(r,t)i  +  y(r,t)j  +  z(r,t)k  .  (9) 

The  functions  x(r,t),  y(r,t)  and  z(r,t)  specify  this  co¬ 
ordinate  system  and  thus  are  sufficient  to  calculate  any 
geometric  quantity  needed.  In  particular,  knowledge  of  equa- 
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tion  (9)  allows  one  to  calculate  the  unit  vectors  ef  and  efc 
in  the  directions  of  r  and  t  respectively.  Then  the  surface 
potential  flow  velocity  can  be  written  as 

V  =  Urer  +  Ufce  .  (10) 

One  should  be  careful  to  note  that  the  components  Ur  and  Ufc 
cannot  be  obtained  from  V  described  in  an  orthogonal  frame, 
such  as  equation  (2),  by  simple  projections.  Suppose  that 
is  described  in  some  other  reference  frame  and  the 
components  and  Ut  are  needed.  Then  the  proper  components 
in  the  (r,t)  system  are  calculated  by  the  formulas 
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where  n  =  er  x  efc  . 


Let  (u,v,w)  denote  the  boundary  layer  velocity 
components  in  the  (r,t,n)  directions  respectively  then  the 
boundary  layer  equations  are  the  followings 


The  continuity  equation  is 


It +  rt  If +  Ik  +  uCr +  -  0 
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The  r-momentum  equation  is 
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The  t-momentum  equation  is 
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(12c) 


with  boundary  conditions 


u=v=w=0  at  n=0 

and  (12d) 

u  =  Ur,  v  =  Ut  at  n=°° 

These  equations  have  been  written  in  terms  of  a 
rotating  coordinate  system  with  the  rotation  rate  vector  ^ 
whose  magnitude  is  ft.  The  component  of  ft  in  the  direction 
normal  to  the  surface  is  denoted  by  fin.  The  pressure  P  in 
these  equations  is  the  reduced  pressure 


p  =  p 


(13) 


where  R  is  the  perpendicular  distance  from  the  point  on  the 
surface  to  the  axis  of  rotation,  P  is  the  Bernoulli  pressure 


obtain  from 
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P  =  i  v  •  V 


p  - 


and  V  =  urer  +  utet  re£ecced  to  the  rotating  reference 
frame. 

The  coefficients  of  equations  (12)  are  defined  in 
terms  of  the  surface  geometry  (the  metrics  hf,  hfc  and  g)  as 
follows: 
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Although  equations  (14)  appear  formidable  and  complicated, 
these  coefficients  depend  only  on  the  geometry  of  the  surface 
and  the  coordinate  system  used  to  represent  this  geometry. 
At  the  actual  calculation  stage  it  is  not  difficult  to  write 
a  subroutine  that  evaluates  these  terms  for  whatever  geo¬ 
metric  representation  is  available.  This  part  of  the  calcu¬ 
lation  can  be  kept  isolated  from  the  details  of  the  boundary 
layer  calculation  method. 


The  boundary  layer  equations  (13b  and  c)  are  non¬ 
linear.  Thus  in  any  procedure  for  solving  these  equations  an 
iterative  approach  is  required.  This  iterative  approach  can 
be  implemented  at  the  present  analytical  stage  or  later  at 
the  numerical  stage  after  the  equations  (13)  have  been  ap¬ 
proximated  by  finite  differences.  It  seems  a  little  easier 
to  formulate  the  iterative  scheme  at  this  stage. 
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Nonlinear  terms  are  decomposed  by  Newton  s  method 

in  which  it  is  assumed  approximate  values  of  the  dependent 

variables  are  available  because  of  the  iteration  procedure. 

2  3  u 

Thus  typical  terms  such  as  u  ,  u-^and  uv  are  replaced  by  the 
expressions 


2  —  —2 
u  =  2uu  -  u 


u  —  9u  ,  9u  -  9u 
r  9r  9r  3r 


UV  =  UV  +  VU  -  uv  , 


where  the  overbar  denotes  that  the  quantity  is  known  from  the 
prior  iteraton  or  as  an  initially  guessed  quantity.  Note 
that  the  above  expressions  are  identities  if  u=u  etc. 


Furthermore  it  will  be  easier  to  work  with  the 
equations  in  first  derivative  form.  Thus,  the  new  variables 
a  and  b  are  defined  by  the  equations 
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By  substituting  expressions  (15)  and  (16)  in  the 
appropriate  places  in  equations  (12)  and  rearranging  the 


terms,  the  boundary  layer  equations  can  be  written 
compact  form 
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It  is  undesirable  to  deal  directly  with  equations 
(17).  In  order  to  circumvent  some  of  the  numercial  dif¬ 
ficulties  associated  with  boundary  layer  growth  and  special 
situations  of  stagnation  points,  symmetry  plane  or  attachment 
lines,  it  is  desirable  to  transform  some  of  the  dependent  and 
independent  variables.  Suppose  that  the  t-direction  is  the 
main  downstream  flow  direction  and  that  S  is  the  arclength  in 
the  t  direction.  Furthermore  suppose  that  each  of  the  com 
ponents  of  X  is  rescaled  by  some  function  of  r  and  t  only. 

Thus,  consider  the  transformations 

r=r,  t=t,  s=n//s\T  (18) 

X  =  D(r,t,/v)*Y  (19) 

where  D  is  a  diagonal  matrix. 


Then  equation  (17)  becomes 
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THE  NUMERICAL  METHOD 


A  three  dimensional  rectangular  grid  is  erected  in 
the  (r,t,^)  coordinate  system.  One  cell  of  this  grid  system 
is  shown  in  Figure  2.  In  this  grid  the  r-coordinate  is 
partitioned  as  {r^},  the  t-coordinate  as  { t ^ >  and  the  C 
coordinate  as  The  value  of  any  function  F  at  a  grid 

point  is  denoted  by  F.  .  .  .  It  is  assumed  that  the  main  flow 

1  /  J  f  * 

direction  is  towards  increasing  values  of  r  and  t  so  that  the 

3  —  3 

convective  operator  A*g^:+  *-n  equation  (20)  is  marched 

•f 

foreward  in  r  and  t.  This  means  that  if  values  of  Y  are 
known  at  stations  (i,j,k),  (i,j+l,k),  (i+l,j,k)  for  all 
values  of  k  then  the  values  of  Y  can  be  calculated  at  all  k 
of  station  (i+l,j+l,k).  The  finite  difference  equations  are 
then  obtained  simply  by  a  three-dimensional  trapezoidal  rule 
integration  over  the  box  with  corners  (i,j,k),  (i,j+l,k), 

(i+1, j,k),  (i+l,j+l,k),  (i,j,k+l),  (i,j+l,k+l),  (i+l,j,k+l) 

and  ( i+1, j+l,k+l) .  This  is  called  the  box  method. 

A  simple  way  to  view  this  method  is  to  consider 
that  the  differential  equation  (20)  is  evaluated  at  the 
center  of  the  box.  Thus,  the  value  of  quantity  at  the  center 
of  the  box  is  denoted  by  an  asterisks  subscript,  such  as  F*. 
Thus  the  value  of  Y  at  the  center  of  the  box  is  simply 
obtained  by  averaging  the  values  over  the  eight  corners  of 
the  box.  Similarly,  the  value  of  a  derivative  at  the  center 


of  the  box  is  obtained  by  averaging  the  two-point  differences 
formed  at  four  appropriate  edges  of  the  box. 


In  order  to  simplify  the  notation  the  following 
scheme  will  be  used.  The  triple  index  (i, j,k)  is  replaced  by 
a  double  index  (l,k)  where  1  =  1,2,3,  or  4  corresponding  to 
the  four  points  (i+l,j+l),  (i,j+l),  (i,j)  and  (i+l,j).  Thus, 
by  looking  down  onto  the  (r,t)  plane  a  typical  face  of  the 
box  at  level  k  is  shown  below. 


(i, j+1)  2 


(i,j)  3 


11  (i+1,  j+1) 


etc. 


14  (i+1, j) 


The  numbering  1  increases  by  proceeding  counter-clockwise 
around  the  box. 


Thus,  the  various  finite  difference  formulas  can  be 
written  as  follows; 


**  =  [('$l,k+^l,k+l)  +  (^2,k+?2,k+l+t3,k 


+  *3,k+l+Y4,k+^4,k+l)]/8 
"  <*l,k+*l,k+l)/8  +  Y*a/8 


(22a) 
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We  split  off  the  values  of  Y*_  at  the  points  1  ■  2,3,  and  4 
and  lump  them  into  the  single  term  Y*a  because  these  are 
known  values.  This  procedure  is  followed  for  the  evaluation 
of  all  the  terms  in  the  equation  (20).  Hence: 


3r|*  4Ar.  [(Yl,k"Y2,k}  +  <*l,k+l"Y2,k+l) 


+  (Y4,k-Y3,k>  +  (*4,k+r*3,k+l)] 


4Ar±  (Yl,k-Y2,k}  +  4Ar..  Y*r 


(22b) 
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The  matrices  A,  B  and  C  and  the  vector  T  are  evaluated  at  the 
center,  *  ,  of  the  box  simply  by  evaluating  each  element  of 
the  matrix  at  the  center  using  the  formulas  (22).  These 
matrices  are  then  denoted  by  the  asterisk  subscript. 


The  result  is  that  the  finite  difference  equation 
corresponding  to  equation  (20)  is  the  followings 
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1 ,  k+1 
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with  6  =4 . 


Thus,  the  system  of  algebraic  equations  (23)  for  k-1,2,..., 
Kt-1,  where  KT  is  the  total  number  of  grid  pints  in  the  C - 
direction,  together  with  the  boundary  conditions  can  be 
solved  for  the  K„  unknown  vectors  Y..  . 
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at  I*  4 At j  [(Yl,k  Y4,k}  +  (Yl,k+l_Y4,k+l) 

+  (Y3,k"Y2,k}  +  (Y3,k+rY3,k}] 
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( 25d ) 


By  substituting  the  formulas  (25)  into  equation 
(20)  and  evaluating  the  matrices  A,  B  and  C  and  the  vector  ^ 
at  the  point  *  one  obtains  equation  (23)  with  terms  (24)  and 
6*2  and  with  the  modification  that  Atj  is  replaced  by  2At^  in 
all  the  terms  of  (24). 
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At  any  level  k  at  which  the  finite  difference 
equation  (23)  is  evaluated  either  the  box  or  the  zig-zag  box 

is  set  up  depending  on  the  sign  of  v.  The  resulting  alge- 
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braic  equations  for  solving  the  unknowns  k  =  Yi+1  j+i  k 
for  k=l,...,KT  will  satisfy  the  zone  of  influence  and  depen¬ 
dence  criteria. 

It  should  also  be  noted  that  this  whole  method  can 
be  used  with  only  a  trivial  change  in  indexing  to  march  in 
the  negative  t  direction  for  decreasing  values  of  j  if  the 
main  stream  velocity  is  in  the  +r  and  -t  direction.  In  this 
case  the  indices  1=3,4  simply  change  to  (i,j+2)  and  (i+l,j+2) 
respectively  for  the  box  method  and  (i,j)  and  (i+l,j+2)  for 
the  zig-zag  box  method. 

Finally,  it  should  be  mentioned  that  the  box  method 
can  be  used  as  a  set  up  for  calculating  the  symmetry  plane 
boundary  layers  along  either  r=0  or  t=0  by  simply  setting 
points  1  and  2  the  same  and  points  3  and  4  the  same  (t=0 
symmetry  plane)  or  points  1  and  4  the  same  and  points  2  and  3 
the  same  (r=0  symmetry  plane). 
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STARTING  A  THREE  DIMENSIONAL  BOUNDARY  LAYER  CALCULATION 


All  boundary  layers  begin  at  a  nodal  point  of 
attachment  and  radiate  outwards  from  this  point  like  a  polar 
coordinate  system.  The  flow  at  the  stagnation  point  can  be 
calculated  in  isolation  of  the  rest  of  the  flow  as  shown  in 
Part  I  of  this  report.  Indeed/  nodal  stagnation  point  flow 
serves  as  the  initial  conditions  of  the  rest  of  the  boundary 
layer. 


Because  the  boundary  layer  radiates  outwards  from  a 
nodal  stagnation  point,  it  is  natural  to  calculate  the  full 
boundary  layer  in  a  polar  surface  coordinate  system  whose 
pole  coincides  with  the  stagnation  point.  In  most  cases,  in 
particular  the  present  case,  there  is  at  least  one  plane  of 
symmetry  passing  through  this  stagnation  point  and  it  is 
convenient  to  measure  the  azimuthal  angle  8  from  this  plane 
of  symmetry.  The  stagnation  point  boundary  layer  calculation 
produces  two  components  of  velocity.  One  component  radially 
outwards  in  the  direction  of  the  symmetry  plane  and  the  other 
component  radially  outwards  in  the  direction  perpendicular  to 
the  symmetry  plane.  Prom  these  two  components  the  velocity 
in  any  radial  direction  can  be  obtained  by  simple  rotation. 


m 


s 


The  next  step  in  the  calculation  is  to  advance  the  solution 
along  the  symmetry  plane  9=0  (say  the  upper  side)  or  0=tt  (the 
lower  side).  These  calculations  can  also  be  made  indepen¬ 
dently  of  the  rest  of  the  boundary  layer.  Once  one  these 
calculations  have  been  made  the  boundary  layer  is  known  on  at 
least  two  boundaries  of  the  strip.  These  boundaries  are  r=0, 
6e  [0,  tt]  ;  r>0,  9=0  or  r>0,  0=tr.  Thus,  the  rest  of  the 
boundary  layer  can  be  calculated  by  advancing  a  step  Ar  from 
the  last  radial  position  and  marching,  or  sweeping  in  9  for 
constant  r,  from  one  boundary,  say  9=0,  to  the  other,  say  9=ir 
(or  vice-versa).  This  scheme  is  illustrated  below. 


0 


In  this  schematic  the  dots  are  points  at  which  the  boundary 
layer  is  known.  The  arrow  indicates  the  direction  of  com¬ 
putational  advance  for  fixed  values  of  r  and  the  point  x  is 
the  current  position  at  which  the  boundary  layer  is  to  be 
calculated.  Note  that  in  the  box  scheme,  indicated  by  the 
dotted  box,  the  three  known  points  are  used  to  calculate  the 
unknown  point.  There  is  just  the  right  amount  of  data  to 
make  the  system  determinate. 
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If  the  calculation  method  cannot  be  arranged  so 
that  rectangular  groups  of  four  points  contain  only  one 
unknown  then  the  box  method/  or  the  zig-zag  box,  by  itself 
cannot  be  applied  because  there  would  be  more  unknowns  than 
equations . 

Consider  then  the  formulation  of  the  boundary  layer 
problem  in  a  locally  rectangular  coordinate  system  near  the 
front  stagnation  point.  Suppose  this  system  is  called  (x,y) 
and  that  x=0  (y  coordinate)  is  the  plane  of  symmetry.  Then 
calculation  of  the  stagnation  point  and  plane  of  symmetry 
yields  the  situation  depicted  below. 


By  advancing  Ax  to  the  next  station  it  can  be  seen  that  at 
the  first  step  the  box  will  have  two  points  of  unknowns  and 
is  thus  indeterminate.  Thus,  a  special  procedure  needs  to  be 
used  to  get  the  first  starting  point  when  advancing  to 
station  x+ax.  It  will  be  shown  in  the  next  section  that  for 
wing-like  bodies  it  is  necessary  to  use  a  locally  rectangular 
system  rather  a  polar  system  to  integrate  the  three  dimen¬ 
sional  boundary  layer  equations. 
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CALCULATION  OF  THE.  FLOW  | 

OVER  AN  OBLATE  SPHEROID-POLAR  COORDINATE  I 
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4.1  COORDINATE  SYSTEMS 

The  most  convenient  coordinate  system  for  calcula¬ 
ting  thee  dimensional  boundary  layers  is  a  polar  system  with 
the  pole  at  the  front  stagnation  point.  For  the  oblate 
spheroid  a  polar  system  with  poles  at  z=+l  is  the  following: 

R  =  'J  1-z2  (  t  cos  0i  +  sin0j)  +  zk  (26) 

In  this  coordinate  system  the  curves  of  z=constant  are 
ellipses  in  planes  of  constant  z  and  the  curves  of  B=constant 
are  ellipses  in  meridianal  planes  through  the  z-axis.  This 
coordinate  system  is  nonor thogonal .  However,  for  a  sphere, 
x  =  l,  the  system  is  orthogonal.  The  curves  6  =  0  and  tt  are  the 
intersections  of  the  body  with  the  x,z  plane.  The  curves 
0=+tt/2  are  the  intersections  of  the  body  with  the  y,z  plane. 
The  curves  0=0  and  ^  are  symmetry  lines  of  attachment  in 
general.  The  curves  8=+t>/2  are  symmetry  lines  of  attachment 


for  zero  angle  of  attack. 


The  metric  coefficients  and  unit  vectors  associated 


with  the  ( 
quantities 


where  h  . 
z 


5/  0)  coordinates  system  are  easily  computed.  These 
are  the  following: 

e  =  (-tzcos0i  -  zsinOj  +  Jl-z2k)  /g 
z  z 

eQ  =  (-TsinGi  +  cos0j)/gQ 


-S- 


2  2 

z  A  cos0 


90 


[~2  71  27 

yji  +  A  cos  0 


=  gz4/Cz2 


=  A- 


2 

z  ge 


2 

g  =  -zA  sin0cos0 


.  and  g  are  the  metric  coefficients. 
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The  potential  flow  in  terms  of  these  coordinates 
are  easily  obtained  by  noting  that  from  equations  (1)  and 
(26)  the  relation  between  the  (y,u>)  and  (z,0)  coordinates  is 


y 


cos9 


1 


cosco 


zA/sin  0  + 


2  2o 

z  cos  0 


(27) 


It  must  be  noted  that  the  velocity  components  from  the  (y,w) 
system  must  be  transformed  to  the  (z,e)  system  by  equations 
(11) . 


When  the  spheroid  is  at  an  angle  of  attack  the 
front  stagnation  point  moves  off  the  pole  (y=0,  w=0)  *  (z=-l, 
0=0)  to  the  point  (y=yQ,  w=0)  =  (z=zQ,  0=0).  The  value  yQ  is 
given  by 


(28) 


from  Eq.  (3).  It  is  possible  to  form  a  locally  (distorted) 
polar  coordinate  system  about  the  stagnation  point.  This  can 
be  done  by  starting  with  a  locally  rectangular  system  about 
the  front  stagnation  point.  The  original  elliptic  co¬ 


ordinates  (  y,  oj)  serve  this  purpose. 


The  following  embedded  circle  transformation  with 


pole  at  (yQ,0)  can  be  used  at  the  front  stagnation  point. 


y  =  yQ  +  (n  -1)  yQr  (r+cos4>)  A  ^ 


sinw  =  (n  -1)  yQr  sin<j>A 


(29) 


2  -1 

where  A  =  (l+2rcos<|>+r  ) 


and  n>l.  The  value  of  n=2  is  appropriate  here.  A  schematic 
diagram  of  this  coordinate  system  is  shown  in  Figure  3.  Note 
that  for  r<<l  the  system  reduces  to  the  simple  polar  form 


y  =  yQ  +  yQ(n  -l)rcos<j> 


sinw  =  yQ(n  -l)rcos<t> 


(30) 


Also,  for  r  -  ^  equations  (29)  reduce  to  a  circle  in  the 
(y,w)  plane  given  by 


y  =  nyQCOsri 


sinw  =  nyQSinn 


where  tan^n  =  ^rrtanT^ 


(31) 


This  coordinate  system  may  not  be  used  for  values  of  r>l/n. 
But  note  that  for  r=l/n  the  outer  boundary  of  this  system  has 


already  encompassed  the  nose  of  the  body  (z=-l).  Then  once 
the  calculation  has  reached  the  radius  r=l/n,  a  switch  can  be 
made  to  the  system  defined  by  equation  (26)  which  has  its 
pole  at  z=-l.  The  boundary  layer  data  in  the  (r  f<p )  system 
only  needs  to  be  interpolated  onto  the  appropriate  (z,9) 
points.  Prom  this  point  on,  there  are  no  problems  due  to 
occurence  of  coordinate  system  poles  on  the  rest  of  the  body. 
The  pole  at  z=l  is  not  a  problem  because  the  boundary  layer 
separates  ahead  of  this  point. 

The  details  of  implementing  these  coordinate  system 
need  not  be  dealt  with  here.  It  is  fairly  obvious,  based  on 
the  discussion  in  Section  2,  what  calculations  need  to  be 
done.  These  algebraic  manipulations  are  straight  foreward 
but  tedious. 

4.2  COMPUTATIONAL  RESULTS  ON  THE  SYMMETRY  PLANES 

As  a  first  test  of  the  computational  method  (after 
successful  basic  checkout  on  a  sphere)  the  boundary  layer  on 
the  horizontal  symmetry  plane  (the  z,y  plane,  x=0)  of  a  25% 
thick  (t=0.25)  oblate  spheroid  was  calculated.  Here  the 
(z,0)  system  of  equation  (26)  was  used.  The  symmetry  plane 
considered  is  thus  the  0=0  curve. 


Figure  4  shows  a  plot  of  the  displacement  thickness 
and  effective  blowing  velocity  on  the  leading  edge  of  this 
oblate  spheroid.  What  is  interesting  to  note  is  that  the 
displacement  thickness  6*  actually  decreases  in  the  neighbor¬ 
hood  of  separation  z=0.233  (near  w  =  l.  15tt/2)  .  Also  the 
effective  blowing  velocity  w*  (a  quantity  of  importance  for 
viscous/inviscid  interaction  calculations)  decreases  towards 
separation.  In  two  dimensional  boundary  layers  these 
quantities  diverge  to  infinity  near  separation.  However,  in 
this  three  dimensional  case  there  is  also  a  cross  flow,  as 
depicted  in  Figure  5  on  a  greatly  expanded  scale.  This  cross 
flow  removes  fluid  from  the  side  of  an  elementary  volume  of 
fluid.  Hence  the  growth  of  the  height  of  this  volume  is 
reduced.  The  cross  flow  shown  in  Figure  5  is  at  the  station 
w=107.4°,  just  slightly  ahead  of  separation.  It  must  be 
noted  here  that  the  separation  that  is  being  referred  to  is 
of  the  component  of  the  flow  in  the  z  direction  around  the 
leading  edge.  The  value  of  w=107.4°  is  just  around  the  tip 
towards  the  trailing  edge  side.  Figure  5  is  a  great  circle 
cut  through  the  x-axis. 

Figure  6  shows  two  cross  flow  velocity  profiles. 
One  at  station  w=101.5°  and  the  other  ato)=107.4°.  The  first 
is  just  slightly  beyond  the  the  w =90°  tip  of  the  disk  where 
cross  flow  reversal  begins.  These  reversed  flow  cross  flow 


profiles  clearly  show  the  concentrated  vortex  that  forms  in 
the  boundary  layer  at  the  tip  of  the  disk.  Because  of  the 
symmetry  of  the  disk  at  zero  angle  of  attack  these  vortices 
are  symmetrically  disposed  over  the  top  and  bottom. 

The  first  question  that  comes  to  mind  when  con¬ 
sidering  nonzero  values  of  angle  of  attack,  a*,  is  whether 
one  of  the  vortices  disappears.  This  is  expected  since  it 
seems  that  angle  of  attack  would  thicken  one  of  these 
vortices  and  thin  the  other.  That  is,  the  pressure  on  the 
lower  (windward)  side  of  the  disk  would  become  more  favorable 
and  over  the  upper  (leeward)  side  the  pressure  gradients  are 
more  unfavorable,  i.e.,  tend  to  promote  separation. 

In  order  to  investigate  these  questions,  the  hori¬ 
zontal  plane  of  symmetry  calculation  must  be  abandoned  and  a 
full  three-dimensional  calculation  must  be  performed.  How¬ 
ever,  in  the  full  three  dimensional  case  of  the  disk  at  a 
nonzero  value  of  the  angle  of  attack  a*,  the  vertical  plane, 
(z,x),  is  a  plane  of  symmetry.  The  full  three-dimensional 
calculation  must  begin  at  this  plane  of  symmetry  and  march 
outwards  towards  the  tip  of  the  disk  at  y=+l.  Thus,  since 
the  vertical  plane  of  symmetry  is  required  for  the  3D 
boundary  layer  calculation,  it  is  instructive  to  consider  it 
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Figure  7  shows  the  distribution  of  the  chordwise 
component  (z-component  for  0=0  and  n  )  of  skin  friction  coef¬ 
ficient  C^z  and  the  cross  flow  component  (in  the  direction  of 
3>0,  but  at  0=0)  C^e  of  skin  friction  coefficient.  The 
windward  side  is  the  bottom  (pressure)  side  of  the  disk  on 
which  nothing  of  great  interest  occurs  until  practically  at 
the  trailing  edge.  Over  the  leeward  (suction)  side  of  the 
disk  it  is  interesting  to  note  how  the  cross  flow  reverses 
direction  z.  Note  that  |cfQ|  is  plotted  here!  However,  the 
streamwise  flow  (z-component)  of  the  flow  is  not  even  close 
to  separating.  Thus,  even  over  the  flat  upper  side  of  the 
disk,  three  dimensionality  of  the  flow  is  very  prominent  and 
there  will  be  embedded  longitudinal  vortices  close  to  the 
longitudinal  centerplane  of  the  disk. 


The  aim  now  is  to  calculate  the  complete  three 
dimensional  flow  field  and  this  work  is  described  in  the  next 
section. 


COMPUTATIONAL  RESULTS  OF  THE  FULL  3D  BOUNDARY  LAYER 


Before  attempting  to  calculate  the  three- 
dimensional  boundary  layer  at  an  angle  of  attack  it  was 
thought  desirable  to  calculate  it  at  zero  angle  of  attack 
first.  Thus  we  consider  the  boundary  layer  on  an  oblate 
spheroid  at  zero  angle  of  attack  in  edgewise  flow.  The  (x,0) 


coordinate  system  of  Eq.  (26)  was  used  since  at  zero  angle  of 
attack  the  front  stagnation  point  lies  right  at  the  pole  of 
this  system. 

When  calculating  the  boundary  layer  on  the  spheroid 
at  a  nonzero  angle  of  attack,  a*,  one  has  to  integrate  out¬ 
wards  from  the  6=0  symmetry  plane,  the  vertical,  (x,z), 
plane. 


It  was  found  that  for  spheroids  of  thickness  t>0.6 
the  boundary  layer  could  be  successfully  calculated  in  this 
way.  The  results  will  not  be  shown  because  they  are  rela¬ 
tively  small  perturbations  of  the  sphere  boundary  layer. 
However,  for  t<0.6  the  calculation  advancing  away  from  the 
vertical  symmetry  plane  becomes  unstable  and  cannot  proceed 
around  the  body. 

Figure  8  illustrates  this  phenomena.  The  z- 
component  skin  friction  coefficient,  CfZ»  is  plotted  versus 
the  azimuthal  coordinate  angle  0.  Note:  the  "z-component* 
does  not  mean  a  component  parallel  to  the  z-axis.  When 
referring  to  the  boundary  layer  z  is  a  parametric  surface 
coordinate  that  runs  roughly  in  the  same  direction  as  the  z- 
axis.  The  points  on  the  surface  corresponding  to  the  z 


parameter  are  in  a  one-to-one  correspondence  with  points  on 


■  Ii, 


the  z-axis  that  have  the  same  value  of  z.  In  Figure  9  the 

small  circles  are  computed  values  of  Cfz  when  the  calculation 

method  advances  from  the  vertical  symmetry  plane  towards  the 

tip  of  the  disk  (0*0  to  0*n/2).  The  solid  curve  is  the  Cfz 

distribution  at  the  same  station  z  (cross-section)  calculated 

by  advancing  from  the  horizontal  symmetry  plane  inwards 

towards  the  vertical  symmetry  plane  of  the  disk.  It  can  be 

seen  that  there  is  a  large  discrepancy  in  these  predictions 

of  Cc  .  The  calculation  that  advanced  from  the  vertical 
f  z 

plane  of  symmetry  broke  down  shortly  after  station  z*0.25. 
The  calculation  that  advanced  from  the  horizontal  plane  of 
symmetry  continued  to  separation  that  started  at  about  z=0.25 
and  6=11/2  (i.e.,  just  beyond  the  maximum  diameter  of  the 

body) . 


The  most  important  aspect  of  this  calculation  is 
that  the  method  of  advancing  the  calculation  azimuthally 
around  the  stagnation  point  starting  from  the  vertical 
symmetry  plane  will  not  work  for  disks  of  thickness  x  less 
than  about  0.6.  The  reason  for  this  is  that  in  a  polar 
coordinate  system,  with  pole  through  the  front  stagnation 
point,  the  boundary  layer  is  advanced  downstream  a  distance 
Ar  and  then  swept  azimuthally  around  the  body  either  from  the 
leeward  side  around  the  edge  and  then  towards  the  windward 
symmetry  plane  or  vice-versa.  In  either  case,  this  azi¬ 
muthal  sweep  will  have  to  advance,  at  some  point,  towards  the 
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edge  of  the  disk.  In  this  part  of  the  sweep,  the  external 
inviscid  flow  component  in  the  azimuthal  direction  points 
against  the  sweeping  direction  of  the  calculation  method. 
Even  the  zig-zag  box  cannot  manage  such  a  strong  reversed 
cross  flow.  The  zig-zag  box  method  can  only  cope  with  weak 
reversed  cross  flow  deep  in  the  boundary  layer.  The  method 
cannot  march  against  the  external  inviscid  flow. 

Some  calculations  were  made  by  advancing  from  the 
vertical  symmetry  plane  around  the  azimuth  for  a  disk  at  an 
angle  of  attack  using  the  coordinates  defined  by  (29).  Again 
for  very  thick  disks,  t>0.6,  the  boundary  layer  could  be 
calculated  successfully,  but  not  for  thin  disks  t<0.6.  We 
are  mainly  interested  here  in  very  thin  disks,  t<0.1,  simula¬ 
ting  propeller  blades.  Thus  a  new  method  for  calculating  the 
3D  boundary  layer  had  to  be  developed. 

If  a  surface  rectangular  coordinate  system  is  used 
in  which  one  set  of  coordinate  lines  run  nearly  longitudinal¬ 
ly  down  the  body  (from  z=-l  to  z=l)  and  the  other  set  of 
lines  run  outwards  from  the  vertical  symmetry  plane  towards 
the  tip,  then  these  coordinates  are  more  nearly  aligned  with 
the  external  flow  over  the  body.  In  terms  of  high  aspect 
ratio  wings,  one  might  say  the  coordinate  lines  are  chordwise 
and  spanwise.  A  calculation  method  based  on  such  a  system  is 


Section  5 


CALCULATION  OP  THE  FLOW 

OVER  AN  OBLATE  SPHEROID-SURFACE  RECTANGULAR  COORDINATES 


5.1  THE  COORDINATE  SYSTEM 

The  coordinate  system  that  is  suitable  for  thin 
bodies  in  nearly  edgewise  flow,  such  as  wings  and  blades,  is 
one  with  coordinates  that  are  roughly  chordwise  and  spanwise. 
Such  surface  coordinates  are  roughly  rectangular  except  for 
the  stretching  required  to  fit  them  over  a  curved  surface.  A 
rectangular  system  does  not  have  any  poles  and  thus  in  the 
neighborhood  of  the  origin  (or  any  other  point)  the  metric 
coefficients  are  approximately  constant.  In  a  polar  co¬ 
ordinate  system  one  metric  coefficient  is  always  a  function 
of  a  coordinate  to  first  order. 

A  suitable  coordinate  system  to  deal  with  the  disk 
boundary  layer  is  the  (r,t)  system  which  gives  the  following 
disk  representation: 


R(r,t)  =  >/l-r2  (xsin  (t-t  )i  +  cos(t-tQ)&)  +  r]  (32) 


where  tQ  corresponds  to  the  value  of  t  of  the  front  stagna¬ 
tion  point  uQ.  By  comparing  equations  (1)  and  (32)/  it  can 
be  seen  that 


y  -  y/l-r^  sin(t-tQ) 


sinw  =  r/vcos^ ( t-t  ) +r  sin  (t-tQ) 


Thus  the  point  r*0/  t=tQ  of  the  front  stagnation  point  cor¬ 
responds  to  o ,  w=0  so  that 


t  =  -sin  1(u  ) 
O  o 


The  coordinate  system  (r,t)  is  nonorthogonal  with 


metric  coefficients 


h  =  >/l-X2r2sin2 (t-t) // 1-r2 
r  o 


hfcv  =  A-X2cos2  (t-tQ) 


g  =  X2rsin (t-t  ) cos (t-tQ) 


The  unit  vectors  e  and  e  are  easily  computed  using  equation 

r  t 


(32)  and  thus  the  potential  flow  velocity 


v  "  ur  er  +  Ut  e  t 
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is  also  easily  computed  with  the  help  of  equations  (3)  and 
(11) . 


It  is  to  be  noted  that  the  surface  coordinate  (r,t) 
has  the  configuration  near  the  front  stagnation  point 
depicted  schematically  below: 

t 


Schematic  of  the  (r,t)  coordinate  system. 


In  this  schematic  0  is  the  front  stagnation  point, 
x  is  the  nose  of  the  disk  on  its  leading  edge  and  the  dashed 
curve  is  the  line  of  attachment  of  the  flow.  The  line  of 
attachment  is  the  curve  on  the  body  which  separates  the 
potential  flow  into  two  regions.  The  flow  above  the  line  of 
attachment  has  a  positive  t-component  of  velocity  that  goes 
over  the  top  (leeward)  side  of  the  body.  The  flow  below  the 


line  of  attachment  has  a  negative  t-component  of  velocity 


that  goes  under  (windward)  side  of  the  body  for  positive 
values  of  the  angle  of  attack  a*. 


THE  CALCULATION  OF  THE  BOUNDARY  LAYER 


The  calculation  of  the  boundary  layer  in  the  (r,t) 
coordinate  system  of  equation  (32)  meets  with  the  difficulty 
of  stepping  the  boundary  layer  from  station  r  to  r+Ar 
described  in  Section  3.  The  procedure  is  basically  to  calcu¬ 
late  the  stagnation  point  boundary  layer  and  the  vertical 
plane  of  symmetry,  r=0,  boundary  layer  for  both  positive  and 
negative  values  of  t.  Once  this  has  been  done,  the  next  step 
is  to  increment  the  value  of  r  to  r+Ar.  At  this  stage  one  is 
faced  with  having  to  compute  the  boundary  layer  at  the  new 
value  r+Ar  without  having  available  any  r+Ar  points  at  which 
the  boundary  layer  is  known.  It  is  clear  from  the  descrip¬ 
tion  of  the  box  and  zig-zag  box  method  of  Section  2.3  that 
these  methods  are  not  applicable  until  at  least  one  point  on 
the  r+Ar  line  (i+1)  station  has  the  boundary  layer  specified. 


A  novel  way  to  obtain  this  required  data  was  first 
suggested  by  Schwamborn  (1981).  In  this  method  the  zig-zag 
box  scheme  is  applied  simultaneously  to  the  first  two  points 
on  the  r+Ar  line  (i+1)  shown  schematically  in  Figure  9.  A 
foreward  zig-zag  box  is  applied  to  the  points  (j+2,i), 
(j+l,i),  (j+l,i+l)  and  (j,i+l)  where  the  point  ( j+1, i+1)  is 


considered  to  be  the  principal  unknown  and  the  point  (i+l,j) 
is  presumed  known  (by  a  guess  or  prior  iteration).  This 
computational  molecules  is  shown  as  the  solid  outline  in 
Figure  9.  A  second,  overlapping,  zig-zag  scheme,  shown  by 
the  dashed  outline  in  Figure  9,  is  used  to  determine  the 
boundary  layer  at  the  point  (i+l,j)  given  the  known  data  at 
points  (i,j),  ( i , j— 1 )  and  (i+l,j+l).  By  simultaneously 
iterating  these  two  zig-zag  schemes  it  seems  that  the 
boundary  layer  at  the  two  points  (i+l,j+l)  and  (i+l,j)  could 
be  determined.  If  this  scheme  converges  then  the  rest  of  the 
points  on  the  line  r+Ar  can  be  calculated  by  the  standard  box 
method . 


Unfortunately,  we  have  not  been  able  to  make  this 
method  work,  although  Schwamborn  claims  that  it  does  work. 
We  believe  the  method  should  converge.  At  this  stage,  our 
inability  to  achieve  success  with  the  method  may  be  due  to 
coding  errors,  though  very  careful  checks,  not  only  by  the 
principal  investigator  but  also  independently  by  another 
scientist  has  not  turned  up  any  coding  errors.  Thus,  the 
method  stands  at  an  impasse  and  no  results  have  been  obtained 
for  the  full  three-dimensional  boundary  layer  calculation  for 
the  oblate  spheroid  at  an  angle  of  attack. 


CONCLUDING  REMARKS 

The  research  effort  to  develop  computational 
methods  for  calculating  propeller  blade  leading  edge  boundary 
layers  is  at  the  following  stage. 

(i)  A  numerical  method  for  dealing  with  the  3D  laminar 
boundary  layer  equations,  with  rotation,  and  in  non- 
orthogonal  surface  coordinates  has  been  developed  and 
programmed. 

(ii)  The  method  has  been  set  up  in  a  a  general  way  so  that 
a  variety  of  coordinate  systems  can  be  handled  easily 
as  demonstrated  in  Sections  4  and  5  of  this  report. 

(iii)  The  method  can  deal  successfully  with  non  flat  bodies 
such  as  ship  bows  and  elongated  bodies  at  zero  or 
nonzero  angles  of  attack,  such  as  ovary  spheroids  or, 
more  practically,  submarines  or  torpedoes.  In  fact, 
the  same  basic  set  of  algorithms  and  codes  were  used 
for  the  direct  numerical  integration  of  the  attachment 
line  boundary  layer  of  Part  I.  This  type  of  calcula¬ 
tion  has  direct  relevance  to  destroyer  bows. 


iv)  A  new  method  for  advancing  from  the  stagnation  point 
along  the  foreward  curve  of  attachment  near  a  leading 
edge  of  a  flattened  body  such  as  a  wing  or  blade  was 
unsuccessful.  The  method  does  not  seem  to  have  any 
flaws  and  at  this  time  it  is  not  certain  whether 
coding  errors  have  prevented  success  in  using  this 
method. 

(v)  Some  alternative  methods  to  achieve  success  in  cal¬ 
culating  the  leading  edge  flow  might  be  to  use  an 
explicit  backward  differencing  with  respect  to  the  t- 
coordinate  or  to  apply  the  characteristic  box  method. 
Resources  of  the  present  contract  have  not  permitted 
us  to  embark  on  efforts  to  develop  these  alternative 
methods.  It  is  believed,  that  a  modest  further  effort 
can  resolve  the  problem  at  the  attachment  line  and 
result  in  a  very  useful,  and  versatile  three- 
dimensional  boundary  layer  code. 
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Figure  1 


Edgewise  flow  past  an  oblate  spheroid  coordinate 
system  schematic. 


Figure  3.  The  embedded  circle  polar  coordinate  system,  s  is 
the  stagnation  point  and  o  is  the  nose  of  the  body 
and  origin  of  the  (y,w)  coordinate  system. 


Figure  4.  Effective  blowing  velocity,  W^,  and  displacement 
thickness  5.  for  t=0.25,  a  =0. 


Skin  friction  coefficients  C£  and  Cca  on  the  vertical 

f  z  to 

(x,z)  symmetry  plane  of  the  disk  of  thickness  t=0.25 
at  an  angle  of  attack  at=0.1. 

Cr-  I - I  Windward  side 

f  z 

C^2  j - 1  Leeward  side 


1 C r Q  I  1 - — j  Leeward  side 


Figure  8.  The  azimuthal  variation  of  the  skin  friction  coef¬ 
ficient  C _  in  the  z-direction  for  a  disk  with 
f  z 

t=0.5,  a*=0  at  station  z=-0 . 25 . 

| - 1  Calculated  by  advancing  from  edge,  0  = 

to  vertical  symmetry  plane  0=0. 


o  o  o  o  Calculated  by  advancing  from  vertical  sym 
metry  plane  0=0  to  the  edge  0=tt/2. 
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